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ABSTRACT 


In this dissertation, an advanced implementation of the direct 
boundary element method applicable to f ree-vibration, periodic (steady- 
state) vibration and linear and nonlinear transient dynamic problems 
involving two and three-dimensional isotropic solids of arbitrary shape is 
presented. Interior, exterior and half-space problems can all be solved by 
the present formulation. 

For the free-vibration analysis, a new real variable bem formulation 
is presented which solves the free-vibration problem in the form of 
algebraic equations (formed from the static kernels) and needs only surface 
discretization. 

In the area of time-domain transient analysis the BEM is well suited 
because it gives an implicit formulation. Although the integral 
formulations are elegant, because of the complexity of the formulation it 
has never been implemented in exact form. In the present work, linear and 
nonlinear time domain transient analysis for three-dimensional solids has 
been implemented in a general and complete manner. The formulation and 
implementation of the nonlinear, transient, dynamic analysis presented here 
is the first ever in the field of boundary element analysis. 

Almost all the existing formulation of BEM in dynamics use the 
constant variation of the variables in space and time which is very 
unrealistic for engineering problems and, in some cases, it leads to 
unacceptably inaccurate results. In the present work, linear and 
quadratic, isoparametric boundary elements are used for discretization of 
geometry and functional variations in space. In addition higher order 
variations in time are used. 


These methods of analysis are applicable to piecewise-homogeneous 
materials, such that not only problems of the layered media and the soil- 
structure interaction can be analyzed but also a large problem can be 
solved by the usual sub-structuring technique. 

The analyses have been incorporated in a versatile, general-purpose 
computer program. Some numerical problems are solved and, through 
comparisons with available analytical and numerical results, the stability 
and high accuracy of these dynamic analyses techniques are established. 


iii 



LIST OF CONTENTS 


Page 

ACKNOWLEDGEMENTS i 

ABSTRACT ii 

NOTATIONS viii 

LIST OF TABLES X 

LIST OF FIGURES xi 

CHAPTER I. INRODUCTION 1 

1.1. The Need for the Present Work 2 

1.2. Relevant Problems of Engineering Analysis and the 

Scope of the Present Work 4 

1.3. Outline of the Dissertation 7 

CHAPTER II. HISTORICAL BACKGROUND 10 

II. 1. Historical Account of Elasto-Dynamics ll 

II. 2. Historical Development of the Boundary Element Method ... 13 

CHAPTER III. REVIEW OF THE EXISTING WORK ON DYNAMIC ANALYSIS 

BY BEM 16 

III. l. Scalar Wave Problems 17 

111. 2. TVo-Dimensional Stress Analysis 18 

111. 2. A. Transient Dynamics 18 

111. 2 . B. Steady-State Dynamics 20 

111. 3. Three-Dimensional Stress Analysis 22 

111. 4. Free-Vibration Analysis 24 

111. 4. A. Determinant Search Method 24 

111.4. B. Domain Integral Transform Method 25 

CHAPTER IV. ADVANCED TWO-DIMENSIONAL STEADY-STATE 

DYNAMIC ANALYSIS 31 

IV. 1. Introduction 32 

IV. 2. Governing Equations 32 

IV. 3. The Boundary- Initial Value Problems of Elastodynamics ... 33 

IV. 4. Boundary Integral Formulation 34 

IV. 5. Numerical Implementation 38 

IV. 5. A. Representation of Geometry and Functions 39 

IV. 5. b. Substructuring Capability 40 

IV. 5. C. Numerical Integration 40 

IV. 5. D. Evaluation of the Diagonal Blocks of F 

Matrix : 42 

IV. 5. E. Diagonal Blocks of F Matirx for 

Problems of Halfspace having Corners and 

Edges 44 

IV. 5. F. Assembly of System Equation 45 

IV. 5. G. Solution of Equations 47 

IV. 5. H. Calculation of Stresses on the Boundary for 

2D Problems 48 


iv 


LIST OF CONTENTS (continued) 

Page 

IV. 6 . Examples of Applications 50 

IV. 6. a. Dynamic Response of a Rigid Strip on an 

Elastic Half space 50 

IV. 6. b. Dynamic Response of a Machine Foundation 

Qnbedded in an Elastic Half space 54 

IV. 6. c. Dynamic Response of a Wall on an Elastic 
Half-space Subjected to a Time Harmonic 
Lateral Pressure Distribution 55 

IV. 7. Concluding Remarks 56 

CHAPTER V. FREE VIBRATION ANALYSIS OF TWQ-DIMESNIONAL PROBLEMS 59 

V. l . Introduction 60 

, V.2. Governing Equation 60 

V.3. Particular Integral 61 

V. 4 . Boundary Element Formulation 64 

V.5. Eigenvalue Extraction 67 

V. 6 . Advantages of the Proposed Method 67 

V.7. Examples of Applications 68 

V.7.a. Comparison with Nardini and Brebbia 68 

V.7.b. Comparison with Finite Element and Beam Theory ... 69 

V.7.c. An Example of a Shear Wall 70 

V.7.d. An Example of an Arch with Square Openings 70 

V. 8. Concluding Remarks 71 

CHAPTER VI. ADVANCED THREE-DIMENSIONAL STEADY-STATE 

DYNAMIC ANALYSIS 76 

VT. 1. Introduction 77 

VI. 2. Boundary Integral Formulation 77 

VI. 3. Numerical Implementation 79 

VI. 3. A. Representation of Geomatry and Field 

Variables 79 

VI. 3. B. Built-in Symmetry and Sub-Structuring 

Capabilities 81 

VI. 3. C. Numerical Integration 82 

VI. 3. D. Calculation of Stresses on the Boundary 

for 3D Problems 85 

VI. 4. Examples of Applications 87 

VI.4.a. Cantilever Subjected to End Shear 87 

VI. 4. b. Cantilever Subjected to Harmonic 

Transverse Load 87 

VI. 4. c. Vertical Compliance of a Rigid Square 

Footing 87 

VI. 5. Concluding Remarks 89 

CHAPTER VII. TRANSIENT DYNAMIC ANALYSIS BY LAPLACE TRANSFORM 91 

VII. 1. Introduction 92 

VII. 2. Laplace Transformed Equations of Elastodynamics 92 

VII. 3. Direct Laplace Transform of Boundary Conditions 93 


v 



LIST OF CONTENTS (continued) 


Page 

VII. 4. Numerical Inversion of Transform Domain 

Solution 94 

VII. 5. Examples of Applications 96 

VII. 5. A. Two-dimensional Applications 97 

VII. 5. A. a. Simply Supported Beam Subjected 

to Step Loading 97 

VII.5.A.b. Half-Space under Prescribed Time- 

dependent Stress Distribution 97 

VII. 5. A. c. Semi-Infinite Beam Subjected to a 

Suddenly Applied Bending Moment 99 

VII. 5. B. Three-dimensional Applications 100 

VII. 5. B. a. Cantilever Beam Subjected to Time- 

harmonic Axial Tension 100 

VII.S.B.b. Spherical Cavity in Infinite Space ... 100 

Vll.5.B.b.i. Spherical Cavity under 

Sudden Radial Pressure 101 

VII.5.B.b. ii. Spherical Cavity Engulfed 

by a Pressure Wave 101 

VII. 6. Concluding Remark 102 

CHAPTER VIII. TIME DOMAIN TRANSIENT DYNAMIC ANALYSIS 103 

VIII. l. Introduction 104 

VIII. 2. Transient Boundary Integral Formulation 105 

VIII. 3. Time-Stepping Scheme 106 

VIII. 3. A. Constant Time Interpolation 107 

VIII. 3. B. Linear Time Interpolation 110 

VIII. 4. Seme Aspects of Numerical Implementation 112 

VIII. 5. Numerical Accuracy, Stability and Convergence 

of Solution 114 

VIII. 6. Examples of Applications 115 

VIII. 6. a. Bar Subjected to Transient End Load 116 

VIII. 6. a. i. Square Cross-section 116 

VIII. 6. a. ii. Circular Cross-section 116 

VIII. 6. b. Spherical Cavity 117 

VIII. 6. b.i. Spherical Cavity under Sudden 

Radial Expansion 117 

VIII.6.b.ii. Spherical Cavity Subjected to a 
Triangular Pulse of Radial 

Pressure 118 

VIII. 6. b.iii. Spherical Cavity Subjected to a 
Rectangular Pulse of Radial 

Pressure 118 

VIII.6.b.iv. Spherical Cavity Engulfed by a 

Pressure Wave 119 

VIII. 6. c. Transient Point Load on Half-Space 119 

VIII. 6. d. Square Flexible Footing on Half-Space 120 

VIII. 7. Concluding Remarks 121 


LIST OF CONTENTS (continued) 


Page 

CHAPTER IX. NONLINEAR TRANSIENT DYNAMIC ANALYSIS 122 

IX. l. Introduction 123 

IX. 2. Boundary Integral Formulation for Dynamic 

Plasticity 124 

IX. 3. Constitutive Model 127 

IX. 4. Discretization and Spatial Integration of the 

Volume Integrals 127 

IX. 4. A. Discretization 127 

IX. 4. B. Spatial Integration 129 

IX. 5. Time-Stepping and Iterative Solution Algorithm 131 

IX. 5. A. Time-Stepping 131 

IX. 5. B. Iterative Solution Algorithm for 

Dynamic Plasticity 133 

IX. 6. Example of Application 134 

IX. 6. a. Bar Subjected to a Step End Load 135 

IX. 7. Concluding Remarks 136 

CHAPTER X. GENERAL CONCLUSIONS AND RECOMMENDATIONS FOR 

FUTURE WORK 137 

X. l. General Conclusions 138 

X.2. Recommendations 140 

REFERENCES 143 

FIGURES 155 

APPENDICES 229 

Al. Boundary Kernels for Two-dimensional Steady-State 

Dynamics A-l 

A2 . Boundary Kernels for Three-dimensional Steady-State 

Dynamics A-4 

A3. Interior Stress Kernels for Steady-State Dynamics A- 5 

A4. Boundary Kernels for Transient Dynamics A-6 

A5 . Interior Stress Kernels for Transient Dynamics A-7 

A6. Volume Kernels for Transient Dynamics A-9 

B Propagation of Wavefronts as Surface of Discontinuity B-l 

Cl. Isoparametric Boundary Elements for 2-D Problems C-l 

C2. Isoparametric Boundary Elements for 3-D Problems C-2 

Dl. Analytical Temporal Integration of the Transient 

Dynamic Kernels for Constant Time Interpolation D-i 

D2. Analytical Temporal Integration of the Transient 

Dynamic Kernels for Linear Time Interpolation D-4 


vii 





displacements and tractions 
stresses 

Kronecker's delta function 
initial stress 

global coordinates of the receiver or field point 

global coordinates of the source point 

displacement and traction fundamental singular solutions 

matrices of coefficients multiplying the known and unknown 


{x}.(y} 

N 


field quantitites, respectively 
known and unknown boundary field quantities 
vector containing past dynamic hisotry 
spatial shape functions for boundary elements 
spatial shape functions for volume cells 


viii 



A incremental quantity 

, spatial derivative 

^ summation 

Superscripts 


. time derivative 

- Laplace or Fourier transformed quantity 

or quantity related to interior stress 

u quantity related to interior displacement 

b quantity related to a boundary point 

s quantity related to elasto-static 


ix 



T.TST OF TABLES 

Page 

4.1. Vertical Stiffness of a Rigid Strip 58 

5 . 1 . Time periods of Free-Vibration of a Triangular 

Cantilever Plate 72 

5.2. Time periods of Free-Vibration of a Square 

Cantilever Plate 73 

5.3. Time periods of Free-Vibration of a Shear Wall 74 

5.4. Free-vibration Modes of Full Arch without and 

with Openings 75 

5.5. Free-vibration Modes of the Syimetric Half of 

the Arch withoutand with Openings 75 

6.1. Comparison of Vertical Compliances Obtained by 

Two Different Meshes 90 


x 



LIST OF FIGURES 

Page 

4.1 Two-dimensional boundary elements 156 

4.2 Boundary element discretization of a half-space problem ... 157 

4.3 Discretization of a rigid strip footing on an elastic 

half-space 158 

4.4 Real part of stiffness coefficients for a rigid 

strip footing 159 

4.5 Imaginay part of stiffness coefficients for a 

rigid strip footing 160 

4.6 Real part of contact stress for vertical vibration 

of a rigid strip footing 161 

4.7 Imaginary part of contact stress for vertical vibration 

of a rigid footing 162 

4.8 Real part of contact stress for horizontal vibration 

of a rigid strip footing 163 

4.9 Imaginary part of contact stress for horizontal 

vibration of a rigid strip footing 164 

4.10 Real part of contact stress for rocking of a rigid 

strip footing 165 

4.11 Imaginary part of contact stress for rocking of a 

rigid strip footing 166 

4.12 Discretization of a machine foundation on 

an elastic half-space 167 

4.13 Real part of stiffness coefficients for a 

machine foundation 168 

4.14 Imaginary part of stiffness coefficients for 

a machine foundation 169 

4.15 Real part of stresses for vertical vibration of 

a machine foundation 170 

4.16 Imaginray part of stresses for vertical vibration 

of a machine foundation 171 

4.17 Real part of stresses for rocking of a machine 

foundation 172 

4.18 Imaginary part of stresses for rocking of a machine 

foundation 173 



LIST OF FIGURES (continued) 

Page 

4.19 A wall in an elastic half-space subjected to 

a time harmonic lateral load 174 

4.20 Lateral displacement of a wall in an elastic 

half-space 175 

5.1 First and fourth bending inodes of a cantilever beam 176 

5.2 Convergence of first six BEM eigenvalues of a 

cantilever beam 177 

5.3 Boundary element discretization of the cantilever beam .... 178 

5.4 Discretizations of a shear wall 179 

5.5 Boundary element discretization of a fixed arch with 

openings iso 

6.1 Three-dimensional nonplanar surface patch 181 

6.2 Three-dimensional surface elements 182 

6.3 Infinite elenent 183 

6.4 Typical subdivision patterns for surface elements 184 

6.5 Typical integration process for a quadrilateral element ... 185 

6.6 Cantilever subjected to harmonic end shear 186 

6.7 Cantilever subjected to harmonic patch load 187 

6.8 Boundary element discretization for a square footing 

on half-space 188 

6.9 Vertical compliance for square footing 189 

7.1 Simple-supported beam subjected to step loading 190 

7.2 Dynamic response of simple-supported beam 191 

7.3a Half-space under prescribed time-dependent 

stress distribtuion 192 

7.3b Time history of displacement u, at the 

internal point F 192 

7.4 Discretization of the half-space 193 

xii 


LIST OF FIGURES (continued) 


Page 

7.5 Time history of displacement u 2 at the 

internal point D 194 

7.6 Time history of displacement u 2 at the 

internal point E 195 

7.7 Time history of displacement u 2 at the 

internal point G 196 

7.8 Stress a 22 at the internal point A 197 

7.9 Stress <r 22 at the internal point B 198 

7.10 Stress a 2 2 at the internal point C 199 

7.11 Semi-infinite beam subjected to a suddenly 

applied bending moment 200 

7.12 Transverse displacement along the semi-infinite beam 201 

7.13 Transient analysis of a cantilever subjected to a 

harmonic axial loading 202 

7.14 Boundary element meshes used in the analysis of 

explosion in a spherical cavity 203 

7.15 Fadial displacement of the cavity surface by 

transform algorithm 204 

7.16 Normalized Hoop stress at the cavity surface by 

transform algorithm 205 

8.1 Time marching process 206 

8.2 Normalized radial displacements of the cavity surface 

by using time steps AT = 0.0002 s. 0.0003 s, 0.0004 s .... 207 

8.3 Normalized radial displacements of the cavity surface 

by using time steps AT = 0.0005 s, 0.0006 s, 0.0007 s 208 

8.4 Normalized ratial displacements of the cavity surface 

by using time steps AT = 0.0008 s, 0.0009 s, o.ooi s 209 

8.5 Normalized radial displacements of the cavity surface 

by using all the three meshes 210 

8.6 Longitudinal stress at the midspan of a cantilever 

beam subjected to an end load 211 


xiii 



LIST OF FIGURES (continued) 


Page 

8.7 Normalized axial displacements at the free end 

of the beam 212 

8.8 Surface discretization of a circular bar 213 

8.9 Normalized axial stresses at the midspan of the bar 214 

8.10 Normalized axial displacements at the free end 

of the bar 215 

8.11 Deviatoric stresses at the cavity surface 

for suddenly applied and maintained pressure 216 

8.12 Radial expansion of a cavity by a triangular 

pulse of radial pressure 217 

8.13 Radial expansion of a cavity by a rectancular 

pulse of radial pressure 218 

8.14 Hoop stresses at the cavity surface for a 

cavity engulfed by a pressure wave 219 

8.15 Radial scattered displacements for a cavity 

engulfed by a pressure wave 220 

8.16 Boundary element discretization for a point 

load on half-space 221 

8.17 Normalized horizontal displacement history 222 

8.18 Transient response of a square flexible footing 

under a prescribed vertical stress disribution 223 

8.19 Disturbance propagation from a point as a 

sequence of co-centric spheres 224 

9.1 Three-dimensional volume cell 225 

9.2 Geometrical mapping of a sub-cell onto a unit cube 226 

9.3 Geometrical mapping of a sub-cell 

(excluding spherical segment) onto a unit cube 227 

9.4 Transient elasto-plastic response of a bar 
subjected to suddenly applied and maintained 

end pressure 228 

xiv 


CHAPTER I 
INTRODUCTION 


1 



1.1 THE NEED FOR THE PRESENT TnPRK 


The dynamic analyses of engineering problems involving two and three- 
dimensional solids have been a subject of intense research for the last two 
decades. For these problems, closed-form analytic solutions are extremely 
difficult to obtain except for very simple geometries and boundary 
conditions which hardly exist in practice. Experiments, on the other hand, 
are expensive and difficult to perform. They also involve elaborate 
apparatus in order to reproduce the desired excitations and to scale the 
important parameters correctly. Therefore, resort has to be made to 
numerical methods of solution. 

There are currently two major categories of numerical methods 
available for dynamic analysis of solids; namely, approximate continuum and 
discrete (lumped parameter) models. The most widely used approximate 
continuum method at present is the Finite Element Method (FEM). In 
principle it appears to be a very versatile technique because it can handle 
complex structure geometry, medium inhcmogeneities and complicated material 
behavior in both two and three dimensions. The finite element formulation 
results in a system of equations that may be solved by modal analysis, 
Fourier transform techniques, or step-by-step integration schemes (Ref. 
Zienkiewicz. 1977). However, the major deficiency of the FEM is that an 
infinite or semi-infinite medium has to be modeled by a mesh of finite 
size. This results in undesirable wave reflections from the artificial 
boundaries. This situation is remedied by the use of transmitting 
boundaries (e.g. Kausel et al, 197 5), hybrid techniques (e.g. Tzong et al, 
1981), or infinite elements (e.g. Bettess, 1977). The use of infinite 
element is restricted to homogeneous far fields because it does not permit 
variation in material properties, and hence problems involving layered 
media cannot be solved by using infinite elements. Similarly, a 
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transmitting boundary encompassing all possible cases of waves impinging at 
the ends of a mesh has yet to be devised. Furthermore, the computational 
cost involved in analyzing three-dimensional problems by the FEM is so 
enormous that only a few researchers can afford it. Another continuum 
method is the Finite difference method (FDM). It has been used less 
frequently than the FEM, primarily because of the difficulties associated 
with it in handling complicated geometries and boundary conditions. 

Discrete models are also in use for a certain class of problems (Ref. 
Hadjian et al, 1974). The basic idea behind the discrete model approach is 
the evaluation of the mass, stiffness and damping coefficients that 
essentially represents the medium. With the use of these frequency 
dependent coefficients known as impedance functions, the dynamic analysis 
of the structure is possible. However, exact expressions for impedance 
functions can be obtained for very few cases only and therefore the use of 
discrete models is rather restricted to some simple problems, e.g. some 
foundation problems (Ref. Arnold et al, 1955; Veletsos, 1971). 

In contrast, it is convincingly demonstrated that accurate and 
efficient solutions to dynamic problems can be easily obtained by using the 
Boundary elenent method (Ref. Banerjee and Butterfield, 1981) because the 
radiation condition is automatically (and correctly) satisfied and for 
linear problems only the surface of the problem needs to be discretized. 
Even for problems with material nonlinearity (e.g. soil), in addition to 
the surface discretization, only a small part of the domain where nonlinear 
behavior is expected needs to be discretized. Thus, a tremendous reduction 
in the size of the problem can be achieved. A brief description of the 
Boundary element method (BEM) is provided in Section II. 2 and a complete 
review of the existing work on dynamic analysis by BEM is presented in 
Chapter III. From this review, it can be seen that most of the existing 
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work on dynamic analysis by BEM suffers either from the lack of generality 
or frcm unacceptable level of accuracy. In addition, all of the existing 
work is based on the assumption of linear elastic behavior and most of them 
assume steady-state conditions. However, in the real world of engineering 
problens, steady-state conditions and linear behavior are at best a first 
order approximation. For truly transient processes it is thus mandatory to 
consider time response and nonlinear behavior. 

Because of the reasons discussed above, there is a need for a complete 
and general analysis method for dynamic problems of two and three- 
dimensional solids, particularly for problems related to the semi-infinite 
mediums. 

The work described in this thesis represents a comprehensive attempt 
towards the development of a general numerical methodology for solving two 
and three-dimensional dynamic problems by using BEM. The developed 
methodology is applicable to free-vibration, periodic vibration and linear 
as well as nonlinear transient dynamic analysis of solid bodies of 
arbitrary shape. 

1.2 KET.EVAOT PROBLEMS OF ENGINEERING ANALYSIS AND 

THE SCOPE OF THE PRESENT WORK 

The ability to predict the dynamic response of solid bodies subjected 
to time and space dependent loads and boundary conditions has gained 
considerable importance in all engineering fields such as machine 
foundation design, seismology, non-destructive testing of materials, soil- 
structure interaction analysis, structural dynamics, metal forming by 
explosives, auto-frettage, and aircraft structure design. 

The methodology for dynamic analysis presented in this dissertation 
can be used for solving a number of problems described above. Brief 
descriptions of some of these problems are given below. 
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(i) Machine Foundation Design : The design of a machine foundation 
essentially consists in limiting its motion to amplitudes and frequencies 
which will neither endanger the satisfactory operation of the machine nor 
will they disturb the people working in the immediate vicinity. Therefore, 
for a successful machine foundation design, a careful engineering analysis 
of the foundation response to the dynamic loads from the anticipated 
operation of the machine is desirable. The existing methods for analyzing 
machine foundations can be categorized into two groups: namely, lumped 
parameter approaches and the finite element method. In the lumped 
parameter approach all the motions are assumed to be uncoupled and for 
complicated geometries it is impossible to find impedance functions. On 
the other hand, as discussed earlier the finite element method is unable to 
handle realistic three-dimensional foundation problems because of its 
finite boundaries and computational costs. Therefore, the methodology 
presented here provides a viable tool for analyzing machine foundations 
with complex geometries embedded in layered soils. The multi-region 
capability of the present code will allow the realistic modeling of the 
foundation as well as the soil. It should be noted that the assumption of 
a rigid or flexible foundation is not needed in the present case. Also, 
different combinations of dynamic loading and boundary conditions can be 
easily incorporated. 

(ii) Seismology : In the field of seismology, one is concerned with the 
study of wave propagation in soils. For this purpose, linearized theory of 
elastodynamics are commonly used. Thus, the present work provides a 
general methodology for studying wave propagation in a homogeneous 
half space as well as in layered soils. 
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(iii) Auto-frettage Process ; This process is used in gun-building and in 
the construction of pressure vessels. In this process, walled structures 
such as pipes and spherical and cylindrical shaped containers are 
deliberately subjected to high pressure during their construction. This 
causes plastic deformation and thereby raises the yield strength of the 
material and induces favorable stress distributions. As a result, the 
working loads (i.e. internal pressures) are now carried out by purely 
elastic deformations. In order to achieve an optimum design of a pressure 
vessel by auto-frettage, the auto-frettage process has to be analyzed 
numerically. For this purpose, nonlinear static analysis algorithms are 
generally used. However, a realistic simulation of this problem can only 
be achieved by using a nonlinear dynamic analysis algorithm. The nonlinear 
transient dynamic algorithm presented in this thesis can serve this 
purpose. 

(iv) Structural Dynamics : The problems related to forced and free- 

vibration of structural components such as beams, columns, and shear walls 
can all be analyzed by the proposed methodology. The nonlinear behavior of 
a structure subjected to an arbitrary transient loading can also be 
obtained by using the present method including the cracking and yielding of 
joints. 

(v) Soil-structure Interaction : The safety of structures such as nuclear 

power plants, dams, bridges, schools, hospitals, and utility pipelines 
during an earthquake is of great concern to the designers and the local 
authorities. Thus, to determine the response of these structures during an 
earthquake, a great deal of research has been done and several techniques 
have been developed. Nevertheless, the problem is so complicated that it 
is still a subject of intensive stucfy. 
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The response of structure during an earthquake depends on the 
characteristics of the ground motion, the surrounding soil, and the 
structure itself. For structures founded on soft soils, the foundation 
motion differs from that in the free-field due to the coupling of the soil 
and structure during an earthquake. Thus, soil-structure interaction has 
to be taken into account in analyzing the response of structures founded on 
soft soils. The available soil-structure analysis techniques can be 
categorized in two groups: i.e., the direct method and the substructure 
approach. In the substructuring approach, one of the steps involved is the 
determination of the dynamic stiffness of the foundation as a function of 
the frequency. The steady-state dynamic algorithm of the present work can 
be used to determine the dynamic stiffnesses of two or three-dimensional 
foundations and embedment of the foundation and layering of the soil can 
both be taken into account. As discussed earlier this methodology is a 
better alternative to the finite element method for this type of problem. 

The time-domain, nonlinear, transient algorithm presented in this 
thesis is a strong candidate for realistic analysis of soil-structure 
interaction problons because, in addition to embedment and layering, it can 
also take into account the nonlinear behavior of soils. Finally, for 
structures subjected to wind load, the present implementation provides an 
accurate and efficient analysis. 

1.3 OUTLINE OF THE DISSERTATION 

This dissertation presents a complete and general numerical 
implementation of the direct boundary element method applicable to free- 
vibration, periodic vibration and linear and nonlinear transient dynamic 
problems involving two and three-dimensional isotropic piecewise 
homogeneous solids of arbitrary shape. 
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A. 


The early history of elastodynamics is presented in Chapter II. Also 
presented is a brief introduction to the boundary element method, its 
historical background and recent developments. 

A literature review of the existing work on dynamic analysis by 
boundary element method is presented in Chapter III. In this chapter, for 
completeness, work on scalar wave problems is also reviewed although it is 
not related to the present work because in elastodynamics waves are 
considered to be vectors not scalars. 

In Chapter IV, an advanced implementation of the direct boundary 
element method for two-dimensional problems of periodic vibrations is 
introduced. The governing equations of elastodynamics are presented 
followed by the boundary integral formulation in transformed domain. 
Subsequently, numerical implementation is introduced which includes 
discussions on the use of isoparametric elements, advanced numerical 
integration techniques, and an efficient solution algorithm. Some 
numerical problems are solved and the results are compared with available 
analytical and numerical results. 

A new real-variable BEM formulation for free-vibration analysis and 
its numerical implementation for two-dimensional problems are presented in 
Chapter V. This method solves the free-vibration problem in the form of 
algebraic equations and needs only surface discretization. First, the 
formulation of the problem is introduced and then some simple problems are 
solved and compared with available results to demonstrate the accuracy of 
this new method. 

In Chapter VI, an advanced implementation of the BEM applicable to 
steady-state dynamic problems of three-dimensional solids is presented. 
The governing equations and boundary integral formulation are the same as 
those introduced in Chapter IV. The numerical implementation for three- 
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dimensional problems is discussed first. Additional features like built-in 
symmetry and sliding at interfaces are also introduced. Finally, a few 
numerical problems are solved and are compared with the available results. 

The Laplace-transform-domain, transient, dynamic algorithm applicable 
to two and three-dimensional solids is introduced in Chapter VII. The 
basic formulation and the inverse transformation techniques are discussed 
first followed by a number of example problems which demonstrates the 
stability and accuracy of this algorithm. 

In Chapter VIII, the boundary element formulation for time domain 
transient elastodynami'cs and its numerical implementation for three- 
dimensional solids is presented for the first time in a general and 
complete manner. Higher order shape functions are used for approximating 
the variation of field quantities in space as well as in time. The 
unconditional stability and accuracy of this algorithm is demonstrated by 
solving a number of problems and comparing the results against available 
analytical solutions. 

Chapter IX presents for the first time in the history of boundary- 
element analysis a direct boundary-element formulation for nonlinear 
transient dynamic analysis of solids and its numerical implementation for 
three-dimensional problems. The formulation is discussed first followed by 
discussions on constitutive model, volume integration, time stepping and 
iterative solution algorithm. Subsequently, a few numerical problems are 
solved and results are presented. 

Finally, conclusions and recommendation for future research are set 
forth in Chapter X. 
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HISTORICAL BACKGROUND 
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II. 1 HISTORICAL ACCOUNT OF ELASTO-DYNAMICS 


The study of wave propagation in elastic solids has a long and 
distinguished history. Until the middle of the nineteenth century light 
was thought to be the propagation of a disturbance in an elastic ether. 
This view was espoused by such great mathematicians as Cauchy and Poisson 
and to a large extent motivated them to develop what is now generally known 
as the theory of elasticity. The solution of the scalar wave equation as a 
potential was first achieved by Poisson (1829). In 1852, Lame added the 
vector potential appropriate to the solenoidal displacement component to 
the Poisson's general solution. Thus, through the efforts of Poisson and 
Lame it was shown that the general elastodynamic displacement field can be 
represented as the sum of the gradient of a scalar potential and the curl 
of a vector potential, each satisfying a wave equation (i.e. longitudinal 
and transverse wave equations). Clebsch (1863), Somig liana (1892), Tedone 
(1897), and Duhem (1898) provided the proof for the completeness of LamSs 
solution; and in 1885 Neumann gave the proof of the uniqueness for the 
solutions of the three fundamental boundary initial value problems for 
finite elastic medium (recently, the proof of the uniqueness is extended to 
infinite medium by Wheeler and Sternberg, 1968). Later, Poisson's solution 
was presented in a more general form by Kirchoff (18 83). This problem of 
scalar wave was further studied as a problem with retarded potentials by 
Love (1904). 

Investigation of elastic wave motion due to body forces was first 
carried out by Stokes (1849) and later by Love (1904). In 1887, Rayleigh 
made the very important discovery of his now well known surface wave. In 
1904 Lamb was the first to study the propagation of a pulse in an elastic 
half-space. He derived his solutions through Fourier synthesis of the 
steady-state propagation solutions. The ingenious technique of Cagniard 
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for solving transient wave problems came along in 1939. He developed the 
technique of solving the problen in the Laplace transform domain and then 
obtained the solution by inverse Laplace transform. This technique is the 
basis for much of the modern work in transient elastodynamics. 

The classical works on elastodynamics are collected and presented with 
the recent analytical developments in a number of books, such as Achenbach 
(1973), Eringen and Suhubi (1975), and Miklowitz (1980). 

During the early I960s > some pioneering work using an integral equation 
formulation was done for acoustic problems by Friedman and Shaw (1962), 
Banaugh and Goldsmith (1963a), Papadopoulis (1963) and others. Kupradze 
(1963) also has done a great deal of work in the extension of Fredholm 
theory to the formulation of problems ranging from linear, homogeneous, 
isotropic elasto-statics to the vibrations of piecewise homogeneous bodies. 
The general transient problem was attempted by Doyle (1966) who used the 
singular solution for the transformed equations to obtain representations 
for the displacement vector, dilatation, and rotation vector. However, he 
did not attack the general boundary value problem in terms of boundary data 
and did not attempt a solution and inversion to complete the problem. 
Nowacki (1964) also treated the transient problem but his solution method 
required finding a Green's function before attempting the Laplace 
inversion. During the past two decades, Banaugh and Goldsmith (1963a) were 
the first ones to use the boundary integral formulation to solve an 
elastodynamic problem. After that, a number of researchers have used the 
boundary element method for solving elastodynamic problems. A complete 
review of these works is presented in Chapter III. 
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II. 2 HISTORICAL DEVELOPMEM* OF THE BOUNDARY ELEMENT METHOD 

The boundary element method (BEM) has now emerged as a powerful 
numerical technique for solving problems of continuum mechanics. In recent 
years, it has been successfully employed for the solution of a very wide 
range of physical problems such as those of potential flow, elastostatics, 
elastoplasticity, elastodynamics, acoustics etc. The BEM, has a number of 
distinct advantages over the Finite element (FEM) and Finite difference 
(FDM) methods such as; discretization of only the boundary of the domain of 
interest rather than the whole domain (i.e., the dimensionality of the 
problem is reduced by one), ability to solve problems with high stress 
concentrations, accuracy, and the ease of solution in an infinite and semi- 
infinite domain. 

This method essentially consists of transformation of the partial 
differential equation describing the behavior of the field variables inside 
and on the boundary of the domain into an integral equation relating only 
boundary values and then finding out the numerical solution of this 
equation. If the values of field variables inside the domain are required, 
they are calculated afterwards from the known boundary values of the field 
variables. The above described transformation of the partial differential 
equation into an integral equation is achieved through the use of an 
appropriate reciprocal work theorem, the fundamental singular solution of 
the partial differential equation (Green's function) and the divergence 
theorem. The BEM yields a system matrix which is much smaller than that of 
a differential formulation (i.e. FEM or FDM) but, in BEM. the system matrix 
is fully populated for a homogeneous region and block banded when more than 
one region is involved. 

Historically, the first use of integral equations dates back to 1903 
when Fredholm (1903) formulated the boundary value problems of potential 
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theory in the form of integral equations and demonstrated the existence of 
solutions to such equations. Since then they have been studied intensively 
particularly in connection with field theory (e.g. Kellog, 1953; Kupradze, 
1963; Muskhelishvili, 1953; Smirnov, 1964). During the 1950s, a major 
contribution to the formal understanding of integral equations was provided 
by Mikhlin (1957, 1965a, 1965b) who studied the singularities and 
discontinuities of the integrands. Due to the difficulty of finding 
closed-form analytical solutions, all of the classical work has, to a great 
extent, been limited to the investigations of existence and uniqueness of 
solutions of problems of mathematical physics, except for the simplest of 
problems (Ref. Morse and Feshbach, 1953). However, the emergence of high- 
speed computers in late 1960s spurred the development of numerical 
algorithms based on adaptations of these integral formulations to the 
solution of general boundary value problems and the resulting technique 
came to be known as the Boundary Element Method. 

The pioneering works in the field of BEM was done by Shaw and Friedman 
(1963a, b) for scalar wave problems; Banaugh and Goldsmith (I963a,b) for 
elastic wave scattering problems; Hess (I962a,b), Jaswon (1963), and Symm 
(1963) for potential problems; Jaswon and Ponter (1963), and Rizzo (1967) 
for elastostatic problems; Cruse (1967) for transient elastodynamic 
problems; Swedlow and Cruse (1971) for elastoplastic problems; and Banerjee 
and Butterfield (1977) for problems of geomechanics. 

In recent years, advances such as the use of higher-order elements, 
accurate and efficient numerical integration techniques, careful analytical 
treatment of singular integrals and efficient solution algorithms have had 
a major impact on the competitiveness of the BEM in routine linear and 
nonlinear two and three-dimensional static analyses. The contributions of 
Lachat and Watson (1976), Rizzo and Shippy, (1977), Curse and Wilson 
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(1977), Banerjee et al (1979, 1985), Banerjee and Davies (1984), Raveendra 
(1984), Telles (1983, 1981), and Mukherjee (1982) should be mentioned. A 
nuirber of textbooks, such as Banerjee and Butterfield (1981), Brefcbia and 
Walker (1980), Liggett and Liu (1983), Mukherjee (1982), Brebbia, Telles 
and Wrobel (1984), and advanced level monographs, such as Banerjee and 
Butterfield (1979), Banerjee and Shaw (1982), Banerjee and Mukherjee 
(1984), and Banerjee and Watson (1986), provide a full description of the 
recent developments in the Boundary element method. 
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CHAPTER III 

REVIEW OF THE EXISTING WORK ON DYNAMIC ANALYSIS BY BEM 
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III.l SCALAR WAVE PROBLEMS 


The phenomenon of scalar wave propagation is frequently encountered in 
a variety of engineering fields such as acoustics, electromagnetic field 
theory and fluid mechanics. The existence of integral equations for scalar 
wave problems in terms of unknown potential functions dates back to 
Kirchoff (1883). However, the use of boundary integral equations to solve 
the scalar wave problems started in early 1960s with Friedman and Shaw 
(1962) solving the transient acoustic wave scattering problem followed by 
Banaugh and Goldsmith (1963b) solving the steady-state (time harmonic) wave 
scattering problem. Since then a number of researchers have contributed in 
this field. Both transient and steady-state behavior have been analyzed 
for wave scattering as well as radiation problems. A radiation problem is 
one where a given displacement or velocity field is specified on a part of 
the surface. A problem wherein an obstacle with a prescribed boundary 
conditions (usually homogeneous) interacts with some incident wave field 
generated by sources elsewhere is called a scattering problem. It should 
be mentioned that both of the above problems are related to infinite or 
semi-infinite space where boundary element method has no competitor. 

Some comparisons of the BEM against the FDM and the FEM are provided 
by Schenck (1967) for time-harmonic, acoustic scattering and radiation, 
Shaw (1970) for transient and time-harmonic, acoustic scattering and 
radiation, Chertock (1971) and Kleinman and Roach (1974) for acoustic 
problems, and Mittra (1973) for the electromagnetic case. For water wave 
problons, the boundary-integral -equation approach has been used by Garrison 
and Seetharama (1971) and Garrison and Chow (1972), with success. Recent 
works on scalar wave problems include that of Shaw (I97 5a,b), Shippy 
(1975), Meyer et al (1977), Morita (1978), Davis (1976), Groenenboom 
(1983), Mansur and Brebbia (1982), and Misljenovic (1982). 
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It should be noted that the scalar wave problem is much simpler than 
the elastodynamic problem because of the reduced dimensionality of the 
parameters involved in scalar problems and because the analytic 
complexities of the fundamental solutions are also not so severe. 

III. 2 TWO-DIMENSIONAL STRESS ANALYSIS 
(A) Transient Dynamics 

The existing work on two-dimensional linear transient elastodynamic or 
visco-elastodynamic problems can be categorized into the following four 
groups. 

(i) Fourier domain solution : In this approach, the time domain 
response is reconstructed by Fourier synthesis of the steady-state 
solutions obtained by a frequency domain BEM formulation. This approach 
has been used by Banaugh and Goldsmith (1963b), Niwa et al (1975, 1976), 
and Kobayashi et al (1975, 1982). Banaugh and Goldsmith solved a problem 
of elastic wave scattering, Niwa et al and Kobayashi et al solved the 
problem of wave scattering by cavities of arbitrary shape due to the 
passage of travelling waves. Kobayashi and Nishimura (1982) also 
introduced a technique for the problems of fictitious eigenfrequency in 
certain exterior problems. 

(ii) Laplace domain solution : This approach involves solution of the 
problem in the Laplace-transform domain by the BEM followed by a numerical 
inverse transformation to obtain the response in the time domain. Doyle 
(1966) was the first to develop a Laplace domain formulation by BEM, 
but he did not solve any problem while Cruse (1967) presented numerical 
results for the two-dimensional problem of the elastic halfspace under 
transient load in plane strain. Numerical results using this approach have 
been also presented by Cruse and Rizzo (1968) and Manolis and Beskos 
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(1981) . 


(iii) Time domain solution : In this approach, the problem is 
formulated in the time domain by the BEM and solved through a step-by-step 
time integration scheme. The fundamental solution used in this approach is 
a function of time and has time retarding properties. This approach has 
been used by Cole et al (197 8) for the anti-plane strain case (i.e. one- 
dimensional problem), by Niwa et al (1980) for the two-dimensional wave 
scattering problem, by Rice and Sadd (1984) for anti-plane strain wave 
scattering problem, and by Spyrakos (1984) for strip-footing problems. 

(iv) Domain integral transform approach : In this approach the domain 
integral related to the inertia term is transformed into a boundary 
integral by approximating the displacements inside the domain. This 
results in a finite element type matrix differential equation formulation 
which can be solved by using a direct time integration procedure such as 
the Wilson theta method, Houbolt method etc. This approach has been used 
by Brebbia and Nardini (1983) to solve a two-dimensional simple frame. 
This method uses a static Green's function instead of time embedded Green's 
functions and therefore it cannot satisfy the radiation condition nor can 
it reproduce the actual transient response at early times. Because of the 
radiation condition, it cannot be used for semi-infinite problems where the 
BEM has a definite edge over all other numerical methods. 

A comparison of the first three approaches on the basis of their 
accuracy and efficiency was done by Manolis (1983). It should be noted 
that, in the above, some simple two-dimensional or anti-plane strain 
elastodynamic problems were solved such as: (a) the case of an unlined or 
lined circular cylindrical cavity under the passage of longitudinal or 
transverse waves; (b) the cases of square or horseshoe shaped cylindrical 
cavities under longitudinal waves; (c) the case of wave propagation in 
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half-planes, etc. 

Host of the above mentioned work suffers from one or more of the 
following: lack of generality, crude assumption of constant variation of 
the field variables in space and time, inadequate treatment of singular 
integrals, and unacceptable level of accuracy. For example. Cole et al 
found the transient dynamic formulation to be unstable, leading to a 
building up of errors as the time stepping progresses; Rice and Sadd found 
that dominant errors in the method arises from integrating the Green's 
function over the singularity and the time domain formulation when applied 
to time harmonic problems reveals solution error propagation; Spyrakos 
finds his flexible strip results to be affected due to the absence of 
corner and edges in his modeling (this is a consequence which arises due to 
the use of constant elements); and Niwa et al (1976) suggest that use of 
higher approximating techniques for time and space variation of field 
variables may improve the accuracy and stability of their method. All 
these fears has been put to rest in the present work by using a higher 
order interpolation function in time and space, taking care of singular 
integral in an accurate and elegant manner (Ref. Sec. IV.4), using superior 
and sophisticated integration techniques and implementing the BEM 
formulation in a complete and general manner. The time-domain transient 
algorithm developed in this work is unconditionally stable and capable of 
producing accurate results for general three-dimensional problems. 

<B) STEADY-STATE (PERIODIC) DYNAMICS 

Two dimensional steady-state dynamic problems have been solved by 
using the BEM by a number of researchers, such as, Banaugh and Goldsmith 
(1963b) and Niwa et al (1975) obtained the steady-state solution of their 
respective problems before reconstructing the transient response by Fourier 
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synthesis. Recently, Dravinsky (1982a, b) used a indirect BEM formulation 
to study two-dimensional problems of plane wave diffraction by subsurface 
topography, Alarcon and Dominguez (1981) applied the direct BEM to 
determine the dynamic stiffnesses of 2D rigid strip footings, and Kobayashi 
and Nishimura (1983) used the direct BEM to obtain steady-state responses 
of a tunnel and a column in the halfspace subjected to plane waves of 
oblique incidence. Askar et al (1984) presented an interesting, 
approximate, iterative boundary-el ement formulation for steady-state wave 
scattering problems which does not require any matrix inversion. He 
presented the results for the problem of wave scattering by a tunnel in 
half-space. Another interesting study has been done by Nakai et al (1984), 
they introduced viscous dashpots in a two-dimensional analysis to simulate 
energy dissipation in the third direction due to radiation. Lately, 
Estorff and Schmid (1984) has applied the BEM to stucfy the effects of depth 
of the soil layer, embedment of the foundation, and percentage of 
hysteretic soil damping on the dynamic stiffness of a rigid strip in a 
viscoelastic soil. Another work related to rigid strip footing was 
recently presented by Abascal and Dominguez (1984, 1985), where they 
studied the influence of a non-rigid soil base on the compliances 
(flexibility) of a rigid surface footing and response of the rigid surface 
strip footing to incident waves. 

In all the works discussed above, the singularity which arises in the 
traction kernels (fundamental solution) is not taken into account properly 
(Ref. Sec. IV.4), and in all of them except that of Kobayashi and Nishimura 
(1983) it is assumed that the field variables remains constant within an 
element. As pointed out by Kobayashi and Nishimura, it is crucial to use 
higher-order boundary elements for boundary modelling of a steady-state 
dynamic problem so that it is fine enough to be compatible with the wavy 
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nature of the solution. In addition, it should be noted that none of the 
above mentioned algorithm, is capable of solving general two-dimensional 
steady-state elastodynamic or visco-elastodynamic problems because they 
cannot take care of corner and edges which are always present in a real 
engineering problem. To remedy all the above discussed problems, this 
thesis presents a versatile steady-state dynamic algorithm by BEM which is 
capable of solving two-dimensional problems involving complicated 
geometries and boundary conditions. 

III. 3 THREE-DIMENSIONAL STRESS ANALYSIS 

Three-dimensional problens of elastcdynamics were not attempted until 
recently principally because of enormous computing requirements and 
formidable task of numerical implementation. In order to reduce the 
computation and complications involved, simplifications of the BEM 
formulation dictated by the nature of the problem to be solved have been 
developed by a number of workers. 

Dominguez (1978a) simplified the steady-state dynamic kernel functions 
for the special case of periodic surface loading on rectangular 
foundations. He also used another simplified formulation (1976b) to study 
the response of embedded rectangular foundations subjected to travelling 
waves. Karabalis and Beskos (1984) have done similar simplifications to 
the time domain transient boundary integral formulation. Yoshida et al 
(1984) used a simplified BEM formulation for determining the response of a 
square foundation on an elastic halfspace, subjected to periodic loading 
and harmonic waves. Tanaka and Maeda (1984) have developed a Green's 
function for two-layered visco-elastic medium, and using this Green 
function in a simplified BEM formulation they numerically calculated the 
compliances for a hemispherical foundation. More complex problems 
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involving the periodic response of piles and pile groups have been 
attempted by Sen et al (1984, 1985a, 1985b), and Kaynia and Kausel (1982). 
They simplified the boundary integral formulation so that only displacement 
kernels are involved in the formulation. Some authors (Ref. Apsel, 1979; 
Dravinski, 1983) have introduced a potentially unstable method involving an 
'auxiliary boundary' so that singular integration can be avoided. In all 
of the above works, the displacements and tractions are assumed to be 
constant within each element. 

Recently, Rizzo et al (1985) and Kitahara and Nakagawa (1985) have 
implemented the BEM formulation for steady-state elastodynamic problems in 
a general form. Rizzo also implemented a mixed-transform inversion to 
obtain the response in the time domain and a technique for the problsn of 
fictitious eigenf reauency in certain exterior problems with homogeneous 
boundary conditions. Kitahara and Nakagawa have introduced a series 
expansion of the periodic kernels for low frequency range, to obtain a 
stable solution at low frequencies. 

In the present work, the direct boundary element formulations for 
periodic dynamic analysis, transformed domain transient analysis and time- 
domain transient analysis have been implemented for problems involving 
isotropic, piecewise-homogeneous, three-dimensional solids. These 
implementations are general and complete in all respects. In addition, for 
nonlinear transient dynamic analysis of three-dimensional solids, the 
direct boundary element formulation and its numerical implementation are 
presented for the first time. To the best of the author's knowledge, a 
comparable system for steady-state and time dependent analyses by the BEM 
has not yet appeared in the published literature. 
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III. 4. FREE -VIBRATION ANALYSIS 


The existing methods for free-vibration analysis by Boundary element 
method can be classified into the following two categories: 

(A) Determinant search method, and 

(B) Domain integral transform method. 

(A) Determinant search method : 

Most of the existing work on the application of BEM to eigenvalue 
problems falls into this category. This includes the work of Tai and Shaw 
(1974), Vivoli and Filippi (1974), DeMey (1976, 1977), Hutchinson (1978, 
1985), Hutchinson and Wong (1979), and Shaw (1979) for membrane (Helmholtz 
equation) and plate vibrations. Niwa et al (1982) also used this method 
for free-vibration problems of Elasto-dynamics. A review of the existing 
work by this approach can also be found in Shaw (1979), and Hutchinson 
(1984). 

In this method, after the usual discretization and the integration 
process, the boundary integral equation for the eigenvalue problem leads to 
a homogeneous set of simultaneous equations, i.e. 

[A(w) 1 (X) = (0) (3.1) 

where the elements of vector (X) are the unknown boundary conditions at 
each node and the coefficients of matrix [A] are the transcendental 
function of the frequency. These coefficients are complex when calculated 
by using the fundamental solution for the corresponding forced vibration 
problem (e.g. Tai and Shaw, 1974; Niwa et al, 1982), or real when 
calculated by using an arbitrary singular solution (e.g. Hutchinson (1978), 
DeMey (1977)). 

The necessary and sufficient condition for equation (3.1) to have a 
non-trivial solution is 
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D = IA(o)) I = 0 


(3.2) 


The eigenvalues are characteristic roots of this determinant. 
However, in the numerical calculation, the eigenvalue can only be 
determined as parameters which attain local minima of the absolute value of 
the determinant. D . as a function of the frequency, u . This requires 
the formation of equation (3.1) for a large number of trial frequencies, 
which makes this method extremely uneconomical for practical applications. 
Moreover, when the eigenvalues are closely spaced, this method may fail to 
give correct eigenvalues. 

As pointed out by Shaw (1979), this approach also leads to fictitious 
roots when an arbitrary singular solution is used rather than a fundamental 
solution. However, Hutchinson (1985) justifies the use of an arbitrary 
singular solution by stating that one can easily sort out the fictitious 
roots fcy a brief look at the mode shapes. 

(B) Domain Integral Transform Method : 

In this approach, the displacements within the domain are approximated 
by some suitable functions. Due to this approximation, the domain integral 
(related to the displacements within the domain) of the integral equation 
is transformed into boundary integrals by using the divergence theorem. 
Since all the integrals of the integral equation are now related to the 
boundary, after some manipulation, the integral equation is reduced to a 
simple algebraic eigenvalue equation. This method was first proposed fcy 
Nardini and Brebbia (1982). A similar way of achieving volume to surface 
integral conversion has also been outlined recently by Kamiya and Sawaki 
(1985) . 

The main advantage of this method is that the boundary integrals need 
to be computed only once as they are frequency independent rather than 
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frequency dependent (as in the case of determinant search method). 
Moreover, since all of the calculations are in terms of real arithmetic, it 
appears to be economical when compared to the determinant search method. 
The method proposed in this thesis has some superficial similarities with 
this method and, therefore, it is briefly reviewed below. 

The governing differential equation for free- vibration of an isotropic 
homogeneous elastic body can be written as: 

2 

^ + P<»> U^ = 0 (3.3) 

where u^ = components of displacement amplitudes 
= stress tensor components 
o> = natural circular frequency 

p = mass density. 

By using the static Kelvin's point force solution the above differential 
equation can be transformed into an integral representation: 

c ij u i(I) = { G i j(x.Dt i (x)ds - J F^x-IJu^xJds 

s s 

+ pw 2 J" u^(a)G^j (z,jL)dv (3.4) 

V 

where x = field point 
JL = source point 
t^ = traction components = v^n^ 
n k = components of outward normal on the boundary 

= traction kernel corresponding to the displacement kernel 

G il 

c ij = 5 ij “ Pij • where p^ is the ]ump term. 

Equation (3.4) not only contains the unknown displacement u^(x) and 
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the traction t^( x ) on the boundary, but also the unknown displacements 
u^(z.) within the domain appearing in the inertial term. In order to 
formulate the problem in terms of the boundary unknowns only, the 
displacements within the domain u^z) are approximated by using a set of 
unknown coefficients and a class of functions f m (z) (superscript m 

denoting the member of the class), such that 

V 2) = a im fm(z) (3 ' 5) 

where 


f m (z) = c - rCz.fjj,) 


(3.6) 


where c = a suitably chosen constant 

^(z*!^) = distance from the point where the function is 
applied to a point z • 

With this approximation, the domain integral of equation (3.5) becomes, 

J U i (z)G ij ( Z .i)dv = 0im J f m (z)G i:j (x.I)dv (3.7) 

V V 


Now if one can find a displacement field with the corresponding 


stress tensor r 


m 

lik 


such that 


T lik,k 5 li f 


(3.8) 


the volume integral in (3.7) can be transformed into a boundary integral 
via the divergence theorem. Thus equation (3.4) can be expressed as (Ref. 
Nardini, 1982). 


CijUi<JL) 


- J G i j(x.I)t i (x)ds + | F i j(x.i)u i (x)ds 
S S 
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- po 2 {-c.j *“(1) + J G ij (x.i)pf i (x)da - | P i;j (x.i) ^(Xidslo^ 

S S 


(3.9) 


where pJJ^ = = trac±ion vector corresponding to the displacement 


field ’ll ^ , where 


m r i-2\i r ~i 9-iov 

li = L5-4^ C + sod-'')! ^l y i “ 90(1-V) 6 li r 


9-10 V . 3 


(3.10) 


After the usual discretization and integration process, equation (3.10) can 
be written in a matrix form as 

[PHu) - [GHt) = pw 2 ( [G] (p) - [FHH)){a] (3.11) 

The relationship between (u) and {a} can be established using 
equation (3.5), i.e. 

(U) - [Q] (a) (3.12) 

where elements of matrix [Q] are simply the values of the functions 
f m (z.) at the nodal points. 

Since matrix [F] is square and possess an inverse, therefore 

(a] = [Q] -1 {u} (3.13) 

It is important to note that [Ql is a fully populated matrix and 
therefore its inversion is costly for a realistic problem. 

Substituting (3.13) into (3.11), we obtain 

[PHu) - [GHtJ = u> 2 [MHu) (3.14) 

where 

[Ml = pUGJlp} - [F] {«/}) [Q]" 1 (3.15) 
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Equation (3.14) is now written in a submatrix form as follows: 


— 




-i 




— — 


'i 

F 11 F 12 


U 1 


G 11 G 12 


fc l 

2 

M 11 M 12 


U 1 
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> 




= to 




F 21 F 22 


U 2 


G 21 G 22 


fc 2 


M 21 ”22 


U 2 

- 


» i 




. t 




l i 


(3.16) 

where u 1 and u 2 are the displacement vectors related to boundaries s^^ 
and s 2 respectively, and t^ and t 2 are the traction vectors related to 
boundaries Sj and s 2 respectively. 

The homogeneous boundary conditions state that on any part of the 
boundary either u or t is zero. Therefore, assuming u 1 = o and t 2 = 0: 

[f 12 Hu 2 } " [G n ]{t i } = 0,2[M i2 Hu 2 } (3.17) 

[F 22 Hu 2 ) - CG 21 3 (tj) = u 2 [M 22 Hu 2 } (3.18) 

From these two sets of equations, £tj) can be eliminated resulting in: 

[FHu 2 ) = u» 2 [MHu 2 } (3.19) 

Equation (3.19) represents the generalized eigenvalue problem. 

Although the method outlined above (first proposed by Nardini and 
Brebbia) eliminates much of the difficulties of the determinant search 
techniques, it still has a number of deficiencies as a practical problem 
solving tool: 

(1) the form of proposed approximation for the internal displacements via 
equation (3.5) seems to be based on a rather ad hoc basis, 

(2) it is rather difficult to find the displacement tensor and the 

corresponding stress tensor t lik to satisfy equation (3.8) for more 
complex problans such as axi-symmetric and three-dimensional ones or 
those involving inhomogeneity and anisotropy. 
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(3) the matrix algebra involved in the construction of the final system 
equations via (3.13), (3.16-18) restricts the method essentially to 
small test problems. In particular, equation (3.19) cannot be formed 
for a multi-region problem where the interface traction and 
displacements must ranain in the system equations for the algebraic 
eigenvalue problem. 

In addition to the two above discussed methods. Benzine (1980) 
presented a mixed boundary-integral finite-element approach for plate 
vibration problems which also reduces the problem to a standard algebraic 
eigenvalue problem. However, his approach is computationally more 
expensive than the Nardini and Brebbia's (1982) method. 
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IV. 1 INTRODUCTION 


In this chapter an advanced numerical implementation of the boundary 
element formulation for the periodic dynamic analysis of two-dimensional 
problems is described. In this implementation, isoparametric curvilinear 
boundary elements are used. The present analysis is capable of treating 
very large, multizone problems by substructuring and satisfying the 
equilibrium and compatibility conditions at the interfaces. With the help 
of this substructuring capability, problems related to layered media and 
soil-structure interaction can be analyzed. 

In the next few sections, the governing equations of elastodynamics 
are presented followed by a discussion on the boundary element formulation 
of elasto-dynamic problems in the transformed domain. Subsequently, 
materials pertaining to the numerical implementation and the solution 
algorithm are introduced. A number of numerical examples are finally 
presented to demonstrate the accuracy of the present implementation. 


IV. 2 GOVERNING EQUATIONS 

The governing differential equation of linear elastodynamics for 
homogeneous, isotropic, linear elastic bodies is called Navier-Cauchy 
equation, which is expressed as 



' )u i, 




+ c„ u. 


2 j,ii 


b : 


u. 


(4.1) 


where u^(x,T) is the displacement vector and bj is the body force 
vector. Indices i and j corresponds to cartesian coordinates; these 
ranges from 1 to 2 for two-dimensional problem and l to 3 for three- 
dimensional problems. Commas indicate differentiation with respect to 
space, dots indicate differentiation with respect to time T , and repeated 
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indices imply the summation convention. 

The constants c^ and c 2 are the propagation velocities of the 
dilatation (P-wave) and distorsional (S-wave) waves, respectively, and are 
given as 

Cj 2 => (X + 2 {jl )/ p c 2 2 = I l/p (4.2) 

where X and p are Lame constants and p is the mass density. 

In equation (4.1) the displacement u^ is assumed to be twice 
differentiable with respect to space and time, except at possible surfaces 
of discontinuity due to shock wave propagations. The kinematical and 
dynamical conditions related to the propagating surfaces of discontinuity 
are discussed in Appendix B. 

Finally, the constitutive equations for the homogeneous, isotropic, 
linear elastic material are of the form 


T ij = P [( V 


" 2c 2 )u m,m 6 


ID 


+ c- 


' (u i,j 


+ u. 


.i» 


where 


5. . 
ID 


is the stress tensor and 
is the Kronecker delta. 


(4.3) 


IV. 3 THE BOUNDARY- TNTTTAT, VAT, HE PROBLEMS OF ELASTODYNAMICS 

For a well posed problem, the governing differential equations (4.1) 
and constitutive equations (4.3) have to be accompanied by the appropriate 
boundary and initial conditions. Thus, the displacements u^x.T) and 
tractions t^(^.T) must satisfy the boundary conditions 

u i (x,T) = q i (x,T) s fcs u 

t^(x,T) = <Jjij( 2 £,T)n^(x) = p^(x»T) 2 ifeS t for T > 0 (4.4) 
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where rij is the outward unit normal at the surface, 

S u is the part of the surface where displacements are specified, 

S t is the part of the surface where tractions are specified and the 
bonding surface of the body is S = S u + s t , 
and the displacements and velocities satisfy the initial conditions: 

U^x/F-O) « O iQ (x) X£V+S 

u^x.T-O) - U iQ (x> X£V+S (4.5) 

In addition, the displacements and velocities have to satisfy the 
Sammerfeld radiation condition at infinity. 

The proof of the existence and uniqueness of the boundary-initial 
value problems of elastodynamics was first provided by Neumann (1995) for a 
bounded region. Later, it is extended to the infinite domain by Wheeler 
and Sternberg (1968). These proofs are discussed in detail in Miklowitz 
(1980, Secs. 1.11 and 1.12), Eringen and Suhubi (1974, Chapter V), 
Achenbach (1973, Sec. 3.2) and Hudson (1980, Sec. 5.3). 

17.4 BOUNDARY INTBCRAT. FORMULATION 

In many practical applications, it is desirable to predict the dynamic 
response of structures under harmonic excitation. If we assume that enough 
time has elapsed after the initial excitation, the transient part of the 
response will vanish and we will be dealing only with the steady-state 
motion. This problem of steady-state motion can be formulated by taking 
the Fourier or Laplace transform of the equations of motion. 

In steady-state, the excitation and response both are harmonic, 
therefore, the displacement and traction will have the form 
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= u^(£,w)e 1 
t.(x,T) = t.(x,a>)e“ 1<1>T 

where w is the circular frequency, 

u^ is the amplitude of the displacement, 
t^ is the amplitude of the traction, and 
i = 


Substitution of (4.6) into the governing differential equation (4.1) 
and cancellation of the common factor e -1<uT yields the Helmholtz equation 


(Oj 2 - = 2 2>u i,ij * c 2 2 “j,ii* l>" 2u j - 0 


(4.7) 


The time variable is thereby eliminated from the governing 
differential equation and the initial-value-boundary-value problem reduces 
to a boundary value problem only. In equation (4.7) the body force is 
assumed to be zero. 

Similarly, substitution of (4.6) in the constitutive equation (4.3) 
and cancellation of the common factor e -luT yields: 


<r. . = p [ ( C - 2 _ 2 c_ 2 )u 5.. + c 2 (U. . + U . .)] 

lj l 2 m.m ij 2 1,] j,i 


(4.8) 


where is the stress amplitude, and is given by 


a. . = t • n 
ID i D 


(4.9) 


Similarly, application of Laplace transform to the governing equation 
(4.1) under zero initial conditions and zero body force, and to the 
constitutive equation (4.3) yields 
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(4.10) 


<C 1 2 - C 2 2> “i,i 

■ • + C^^U . . . 

L3 2 3,11 

“lj * » [(0 1 2 - 

2c 9 2) u m m 5 . . 
2 in # m 




where the Laplace transform f( X ,s) 
T is defined as 



= o 



+ 



(4.11) 

(4.12) 


of a function f(x,T) with respect to 


L(f ( X ,T)} = f( X .s) = | f( X ,T)e -sT dt (4.13) 

o 

where s is the Laplace transform parameter. 

A comparison of equation (4.7)-(4.9) with (4.10)-(4.12) indicates that 
the steady-state, elastodynamic problem can be solved in the Laplace domain 
if the complex Laplace transform parameter s is replaced by -iu> , o> 
being the circular frequency. It should also be noted that the transformed 
Navier-Cauchy equations are now elliptic, and thus more amenable to 
numerical solutions. 

The boundary integral equation in the Laplace transformed domain can 
be derived by combining the fundamental, point-force solution of equation 
(4.10) with the Graffi's dynamic reciprocal theorem (Graffi, 1947), as 


c ij (JL) u i(£»s) = | [G i j( X ,I,s)t i ( X ,s) - F i j( x ,i,s)u i ( X ,s)]dS( x ) 

S (4.14) 

In the above equation, £. and x are the field points and source points, 
respectively, and the body force and initial conditions are assumed to be 
zero. The fundamental solutions G^j and F— (Ref. Cruse and Rizzo, 
196 8) are the displacements and tractions at X , resulting from a unit 

_ 4 mT ST 

harmonic force of the form e (or e ) at 1 and are listed in 
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Appendix Al. It can be seen that these fundamental solutions have modified 
Bessel functions embedded in them. The asymptotic series expansions of 
these functions for small and large values of argument (i.e. frequencies) 
are also discussed in the Appendix Al. 

The tensor c^ of equation (4.14) can be expressed as: 


c. . - 5. . - a. . 
ID i] P iD 


(4.15) 


where is the discontinuity (or jump) term and it has the following 

characteristics: (i) for a point £ inside the body p^j = 0 , (ii) for a 
point £ exterior to the body p = 5^ , and (iii) for a point £ on the 
surface it is a real function of the geometry of the surface in the 
vicinity of £ . For Liapunov smooth surfaces, p^ = 0.5 5^j . 

Once the boundary solution is obtained, equation (4.14) can also be 
used to find the interior displacements; and the interior stresses can be 
obtained from 


Tj k (£.s) = j [G?j k (x.£»s)t^(x.s) - F?jj < (s,£,s)u^(x,s)]dS(x) 


(4.16) 


The functions and F?j k °f the above equation are listed in 

Appendix A3. 

The stresses at the surface can be calculated by combining the 
constitutive equations, the directional derivatives of the displacement 
vector and the values of field variables in an accurate matrix formulation 
(Ref. Sec. IV.5.H) . Also, the loads and moments on the elements can be 
obtained by numerically integrating the known tractions on the elements. 

The boundary integral formulation can also take account of internal 
viscous dissipation of energy (damping); this can be accomplished by 


37 



replacing the elastic parameters X and u (Lame constants) by their 
complex counterparts x* and n* • 

x* = x(l + 2ip) 

ti* = n(l + 2ip) (4.17) 

leaving Poisson's ratio unaltered. By analogy with single degree-of- 
freedom systems, the damping ratio p is equal to ui\/2\i , where n is 
the coefficient of viscosity for a Kelvin-Voigt model. 

17.5 NUMERICAL IMPLEMENTATION 

The boundary integral equation (4.14) cannot be solved analytically 
and therefore resort must be made to the numerical methods of solution. 
The basic steps involved in a numerical solution process for the boundary 
element formulation are: 

(i) Discretization of the boundary into a series of elements over which 
the geometry and the variation of displacements and tractions are 
approximated by using a suitable set of shape functions. 

(ii) Application of the equation (4.14) in discretized form to each nodal 
point of the boundary and thereby evaluation of the integrals by a 
numerical quadrature scheme. 

(iii) Assembly of a set of linear algebraic equations by imposing the 
boundary conditions specified for the problem. 

(iv) Finally, the system of equation are solved by standard methods to 
obtain the unknown boundary tractions and displacements. 

In the present work, the numerical implementation of the transformed 
boundary element formulation for two-dimensional problems of elastodynamics 
has the following aspects and features: 


38 


(A) Representation of Geometry and Functions 

For the discretization of equation (4.14) the boundary S is 
approximated by using a series of elements whose geometry is defined using 
the quadratic shape functions of intrinsic coordinates proposed by 
Ergatoudis (1968). The boundary elements for two-dimensional problems are 
shown in figure 4.1. On each element the variation of the cartesian 
coordinates x^) are approximated as: 

*i(n) = (4.18) 

where X. are the nodal coordinates of the element, and N are the 

«. let a 

interpolation functions (Ref. Appendix Cl). For a quadratic variation a 
ranges from 1 to 3, and for a linear variation it ranges from 1 to 2. 

Isoparametric shape functions are used to approximate the variation of 
displacements and tractions over each element. In some cases, the full 
quadratic variation of the field quantities is not required so the option 
of using the linear, the quadratic or a mixture of linear and quadratic 
interpolation functions for displacement and traction variation is 
provided. However, the boundary is always- modeled using the quadratic 
shape functions. Using the interpolation functions, the displacement and 
traction at an arbitrary point of a boundary element are expressed in terms 
of nodal values of displacements and tractions by: 

V*) " N a (l t )u ia 

ti ( T,) = N a (n>t ia (4.19) 

where n is the intrinsic coordinate which ranges from 0 to l, and 

u^ q and t^ a are the values of the displacement and traction vectors 
at node a . 
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(B) Substructurina Capability 


In the present implementation, the substructuring capability is 
provided. This is a very useful tool for solving problems related to 
piecewise homogeneous material, layered media and soil-structure 
interaction. This technique actually allows a problem geometry to be 
modelled as an assembly of several generic modeling regions (GMR). The 
GMRs are joined by enforcing appropriate compatibility conditions across 
common boundary elements. 

(C) Numerical Integration 

Taking into account the boundary discretization and function 
representation, the transformed, boundary- integral equation (4.14) can be 
written as: 

Q 

CydJuja.s) - 5 [J [G ij (x(n).i.a)N a (n)t ia ds(x(n)) 

^ S q 

- J P i3 (x(n).i,S)N 0 (n)U ia (x(n))] (4.20) 

S q 

In the above equation, is the length of the qth el orient and Q is 
the total number of elements. In order to express dS(&) in intrinsic 
coordinates, we have 

dS(x> = IJl dii (4.21) 

where 1 J I is the Jacobian which performs the transformation from the 
cartesian coordinate system (x,y) to the elements intrinsic coordinate 
system r\ , and is given by 
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|j(n)l = [ 


(4.22) 


<3N q (h) 

dn 



dN a (n) 

dn 



Therefore, in view of the above, the equation (4.20) can be written as 

Q A t 

c^Du.d.s) = J 5 ha I Gij <x(*>. A.sm o (n)IJldn 

q=l a=l O 
A _ 1 _ 

- 5 u ia J P i j<X(n).i.S)N a (n)ljldn (4.23) 

a=l O 

where A is the number of nodes in an element. 

The global system of boundary element equations is obtained by the 
usual nodal collocation scheme, i.e., by allowing field point I in 
equation (4.23) to coincide sequentially with all the nodal points of the 
boundary. All the boundary integrals involved are calculated numerically. 
Essentially two types of integrals, singular and nonsingular, are involved. 
The integrals are singular if the field point for which the equations being 
constructed lies on the element being integrated. Otherwise, the integrals 
are nonsingular although numerical evaluation is still difficult if the 
field point and the element being integrated are close to each other. 

In both singular and nonsingular cases a Gaussian quadrature scheme is 
used. The basic technique was first developed by Lachat (1975) and is 
discussed in detail by Watson (1979) and Banerjee and Butterfield (1981). 
For the nonsingular case, an approximate error estimate for the integrals 
was developed by Lachat based on the work of Stroud and Secrest (1966). 
This allows the determination of element subdivisions and orders of 
Gaussian integration which will assure roughly uniform precision of 
integrations throughout the integration process. In the present work, this 
automatic choice of integration order and element subdivision has been 
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implemented; where the order of integration points varies from 2 to 12 and 
the number of element subdivisions varies from l to 4. when the field point 
is very close to the elonent being integrated, use of a uniform subdivision 
of the element leads to excessive computing time. Therefore, in order to 
improve efficiency while still retaining accuracy, a graded element 
subdivision is employed. This subelement division grows geometrically away 
from the point closest to the field point on the element being integrated. 

In the case of singular integration, which arises when the field point 
is on the element being integrated, the elementsis divided into 
subelements. The nature of this division depends on the node of 
singularity of the element. This division produces nonsingular behavior in 
all except one of the required integrals. Normal Gaussian rules are used, 
with orders 4 to 8. The integral of the traction kernel times the shape 
function which is l.o at the source point is still singular and cannot be 
numerically evaluated with reasonable efficiency and accuracy. Hence, this 
integral is evaluated indirectly by a scheme discussed in the next section. 

The integration of the surface integrals required for the calculation 
of displacement and stress at interior points are carried out in the same 
manner as that for boundary values (described above) except, in this case, 
all the integrals are nonsingular. 

(D) Evaluation of the Diagonal Blocks of F Matrix 

The diagonal 2x2 block (or 3x3 block for three-dimensional problems) 
of the assembled P matrix contains the tensor c^ as well as the Cauchy 
principal value of the traction kernel integral, i.e. 

5 ij - «« + S Vi as < 4 - 24 > 
s. 
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where: 


is the term which depends only on the geometry at the 
singular node, 

D^j is the diagonal 2x2 (or 3x3 for 3D) blocks of the assembled 
F matrix for the dynamic problem, 

Fj, j is the singular traction kernel for the dynamic problem, 

Nj is the shape function for the singular node, and 

Sj is the length of the singular element. 

Similarly for a static problem: 

% ■ c ij * 1 p !j N i ds ’ (4 - 25> 

s i 

where the variables are the static counterpart of those of equation (4.24). 
From (4.24) and (4.25) we can obtain 


D. . = D?. + f (F- . - F?.)N, 
l] rj J 13 i] l 


dS 


(4.26) 


In the above equation, the diagonal blocks of coefficients of the 
traction ratrix,for a static problem of the same geometry can be obtained 
by using the rigid body motion, i.e. 

°ij - c ij * I 

S, 


Q A 


■ - [ H P lj N a dS P ?jV S ] 


a=2 S, 


q=2 a=l S. 


(4.27) 


In addition, the integral involving the difference (F^j - F^j) is 
nonsingular, therefore, equation (4.26) can be used to obtain D^j . 
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Recently, a somewhat similar approach is used by Rizzo et al (1985) for 
three-dimensional problems. 

In almost all of the past works, the nonsingular integral of equation 
(4.26) has been neglected. This results in inaccuracy, particularly at 
high frequencies. However, for problems related to ground surface the 
above technique is not applicable. Thus, for half space problems a new 
scheme is developed to calculate the diagonal blocks of F matrix. This 
scheme is discussed in the following section. 

(E) Diagonal Blocks of F Matrix for Problems of Half space 

Having Corners and Edges 

The conventional approach of assuming 0.55 „ as the block diagonal 
terms of the F matrix does not hold true for cases where the geometry of 
the problem has corners and edges except for the case where the field 
variables are assumed to be constant within each element. Thus, for higher 
order variation of the field variables, one needs to have a general method 
for calculating the diagonal blocks of F matrix for half space problems. 
In the present work, a new technique to handle the above discussed problem 
in an approximate manner has been developed. To this purpose, this new 
technique uses special types of elements called 'enclosing elements' (Ref. 
figure 4.2). The basic assumption in this technique is that the 
displacements and tractions at the enclosing elements has negligible effect 
on the displacements and tractions at any point on the modeled boundary. 

Using this scheme, the diagonal blocks of F matrix are 

obtained by the summation of nonsingular integrations of the static 
traction kernel over all the boundary elements as well as all the enclosing 
elements, i.e. 
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A Q A 

“ij ' - [ l I F ij N a dS + l H P ij H a aS 

a=2 S 1 q=2 a=l S 

A '"3 

L A 

*5 H P ?j N a aS ] (4 - 28 > 

e=i a=l S o 
e 

where the third summation of the integrals corresponds to the enclosing 
element (L being the total number of enclosing elements). Once D? ^ is 
evaluated, the diagonal blocks related to the dynamic problem can be 

easily found by using equation (4.26). 

In order to show the validity of the above technique, the dynamic 
response of a rigid strip on an elastic halfspace under vertical loading is 
analyzed by using this approach and other two approaches. The real and 
imaginary part of the vertical stiffness for two different frequencies were 
tabulated in table 4.1 obtained by using all the three approaches for 
calculating the block diagonals of F matrix. It can be seen that results 
obtained by using the enclosing element technique compares well with the 
correct results (method 1). However, method 2 which is invariably used by 
the past researchers gives erroneous result at high frequencies (e.g. 
compare the real part of the stiffness at non-dimensional frequency a Q = 
<ob/c 2 = 7.0) 

(F) Assembly of System Equations 

Once the boundary collocation and integrations are completed, we have 
a set of coefficients which function as multipliers of field quantities, 
i.e. (Ref. Banerjee and Butterfield, 1981): 

EG] (t) - [FHuJ - {0} (4.29) 

where: 
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[G] is an unassembled matrix whose coefficients are the values 
obtained by the numerical integration of the product of the 
tensor G^j , the shape functions and the Jacobian. The 
size of matrix CGI is dn x dm ; 

[Fl is an assembled (for nodes) matrix whose coefficients are 
obtained by the numerical integration of the product of the 
tensor F^ , the shape functions and the Jacobian. The 
size of matrix IF] is dn x dn ; 

ft) and Cu } are the transformed traction and displacement 

vectors at the boundary of the problem, with size dm and dn, 
respectively ; 

n is the total number of nodes; 

Q 

m = ^ Ag , where Q is the total number of elements and 
q=l 

Aq is the number of nodes in the qth element; and 

d is the dimensionality of the problem (i.e. for two- 
dimensional problems d = 2) . 

Since some of the field quantities are known from the specified 
boundary conditions, during the assembly of the system equations the 
coefficients related to the known and unknown variables are separated. For 
the case when the boundary conditions are specified in local coordinate 
system, the corresponding coefficients of the matrices [G] and [Fl are 
multiplied by the appropriate local transformation matrix. Finally, 
boundary conditions are imposed including ary required modification to the 
coefficient matrices for bonded or sliding contact between different 
regions (GMRs). The results of all the above operations is a linear system 
of matrix equations of the form: 
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[A] {x} = [BHy} = (b) 


(4.30) 


(u in ) = [A u Hx) + [B u ] Cy } (4.31) 

{<y} = [A®Hx) + [B ct ] {y ) (4.32) 

where 

{y} and (x) are the vectors of known and unknown field 
quantitites, respectively; 

(u in } and { er) are the vectors of displacements and stresses at 
interior points, respectively. 

In any substructured (multi-zone) problem, the matrix [A] in (4.30) 
contains large blocks of zeros because separate GMRs communicate only 
through common surface elements. In order to save both storage space and 
computer time, the matrix [A] is stored in a block basis with zero blocks 
being ignored. Since interior results in any GMR involves only the 
boundary values related to that GMR, the matrices in (4.31) and (4.32) are 
also block diagonal. In addition, for added accuracy the system equations 
are scaled so that all the coefficients of matrix [A] (and CBl) are of the 
same magnitudes (for detail. Ref. Banerjee and Butterfield, 1981). 

(G) Solution-Of Equations 

Since the system equations (4.30) are complex it requires a complex 
solver. In the present work, an out-of-core complex solver is developed 
using softwares from LINBACK (Dongarra et al, 197 9). In this solver in 
order to minimize the time requirements the solution process is carried out 
using block form of the matrix. Thus, this block banded solver operates at 
the submatrix level using software from LINPACK to carry out all operations 
on submatrices. The system matrix is also stored by submatrices on a 
direct access file. The first operation in the solution process is the 
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decomposition of the system matrix using the block form of it. This 
decomposition process is a Gaussian reduction to upper triangular 
(submatrix) form. The row operations required during the decomposition are 
stored in the space originally occupied by the lower triangle of the system 
matrix. Finally, the calculation of the solution vector is carried out by 
using the decomposed form of the system matrix from the direct access file. 


Once the boundary solution is obtained, the stress and strain at any 


point on the boundary can be calculated without airy integration by using 


the procedure outlined as follows. 

Let us assume that we are interested in finding stress and strain at a 
point P , which lies on a boundary element and has intrinsic coordinate 
ii b . Eecalling equations (4.19), we can write: 



where: 

A is the number of nodes in the element, 

N q is the shape functions, and 

u. and t- are the nodal values of u_- and t. . 
la la 11 
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(4.35) 


; U k.l +U l,k. 
*ij = C ijkl ( 2 5 


u. = u. .X, 
i.n i»d 


* 3N 

u. - y r- 2 u. 
i»n 4. 3n la 

0=1 


(4.36) 


(4.37) 


where: 


c ijkl 



is a tensor containing elastic constants, and 
are the directional derivatives. 


Equations (4.3 4), (4.3 5) and (4.36) can be combined to form a matrix 

equation: 
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where n, and n 2 are the unit normal on the boundary at point P ; i.e. 


n = X_ and n_ = -X. 

1 i, J. ,TJ 


Now, the stress and strain at point P can be obtained by inverting 
the matrix of equation (4.3 8) and then multiplying the inverted matrix by 
the right- hand-side vector. For .this purpose, the right hand side vector 
is obtained by using equations (4.33) and (4.37). The procedure described 
above is valid for both plane stress and plane strain problems. However, 
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for plane strain problems the Poisson's ratio v has to be replaced by \> 

V 

l+v • 

17.6 EXAMPLES OF APPLICATIONS 

In order to demonstrate the accuracy and applicability of the present 
implementation, the detailed solution of three numerical examples are 
presented. In the first example, the forced oscillations of a massless 
rigid strip foundation on an elastic halfspace (plane strain) subjected to 
external dynamic forces is analyzed. The purpose of this example is to 
compare the response predicted by the present implementation to that 
available in the literature. The second example is that of a machine 
foundation embedded in an elastic halfspace (plane strain) and subjected to 
external dynamic forces, and the third example is a wall in an elastic 
half-space subjected to a time-harmonic lateral pressure distribution. 
These last two examples are intended to show the applicability of the 
present implementation to real engineering problems. In both examples, 
English units are used with foot (ft.) for length, pound (lbf.) for force, 
and second (s) for time. 

(a) Dynamic Response of a Rigid Strip on an Elastic Halfsrace 

A large number of numerical results have been published for the rigid 
strip with vertical, horizontal and rocking vibrations (Karasudhi et al, 
1968; Luco et al, 1974; Luco and Westmann, 1972; Wickham, 1977; Hryniewicz, 
1981; etc.). However, most of them are limited to a small range of 
frequency parameter and are based on the assumption that one of the contact 
stress components is zero. For the purpose of comparison, a rigid strip 
footing on an elastic halfspace under relaxed boundary conditions is 
analyzed for vertical, horizontal and rocking vibrations. The rigid strip 
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footing and the boundary element mesh used are shown in figure 4.3 (this 
mesh was selected after a convergence study at a high frequency). In all 
cases, a homogeneous soil material with a Poisson's ratio \> = 1/4 is 
considered. The dynamic stiffnesses obtained by the present method are 
compared with that reported by Hryniewicz (1981). He defines the vertical, 
horizontal and rocking stiffness coefficients by the following expressions 
respectively: 


K il + lC il 


P 

miWo 


(4.39) 


K 22 + iC 22 


H 

nuu 0 


(4.40) 


where : 


K 33 + iC 33 


H 

n|ib 2 d 0 


(4.41) 


P, H and M are the amplitudes of vertical force, horizontal force 
and moment, respectively; 

w Q , u Q and <6 Q are amplitudes of vertical displacement, horizontal 
displacement and rotation, respectively; 

K ll' K 22 K 33 are rea ^ P arts t ^ ie stiffness coefficients; 

c ii’ C 22 an ^ C 33 are i ma 9i nar y parts of the stiffness 

coefficients; 

2b is the width of the footing; 

n is the shear modulus of the soil; and 


i =7-1 


The real part of the stiffness coefficients are plotted against non- 
dimensional frequency (a Q = wb/c 2 » where u> is the excitation frequency) 


51 



in figure 4.4. It can be seen that the present results are in good 
agreement with the results of Hryniewicz (1981) for low to medium 
frequencies. However, for higher frequencies the the agreement is not 
good, particularly for rocking stiffness. This difference is due to the 
fact that in the present work quadratic shape functions are used for 
representation of the variation in the boundary variables over each element 
whereas Hryniewicz assumes that the unknown contact stresses are constant 
within each element. This results in stress discontinuities at the 
interface of two elements. Therefore at high frequencies, Hryniewicz's 
method will produce correct results only when the foundation-soil interface 
is divided into a very large number of elements. Figure 4.5 shows the plot 
of imaginary part of the stiffness coefficients against the nondimens ional 
frequency a Q . A good agreement between the present results and the 
results due to Hryniewicz can be seen. Real and imaginary parts of 
vertical stiffness for a bonded rigid strip are also plotted in figures 4.4 
and 4.5, respectively. The imaginary part is identical to that of a 
frictionless rigid strip. 

Dynamic contact-stress distributions at the interface between the 
rigid strip and the halfspace are also presented. For the purpose of 
plotting, the contact stresses are defined as follows: 

For vertical vibration : 


°zz 




+ * CT zz )e 


iwT 


xz 


7TUW „ 

-j-5 (7 1 

b xz 


* 1 


ib/I 


(4.42) 
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For horizontal vibration : 


niIu o ,-r . -i . ion; 

°zz ‘’zz + 1 “zz )e 


a xz 


*^0 -R . -I ia/r 
“ (<J xz + 1 a xz )e 


(4.43) 


For rocking : 


zz 


, R , . I . iwT 
- ntia5 o (ff zz + 1 °zz )e 


a = - m id) (e? + i a 1 

xz ^ o xz xz 


(4.44) 


where superscripts R and I represent real and imaginary parts, 
respectively. 

The real and imaginary parts of the contact stress distribution for 
vertical vibration are plotted in figures 4.6 and 4.7, respectively. 
Because of the singularity at the edge, the contact stresses on the element 
close to the edge are obtained in an average sense (by taking the average 
of nodal values) and are indicated by dashed lines. From the figures, it 
can be seen that the contact stresses are quite sensitive to variations in 
the frequency parameter a Q . As frequency increases, the imaginary part 
of the contact stress distribution increases and the singularities at the 
edge gets sharper for real and imaginary parts. Figures 4.8, 4.9, 4.10 and 
4.11 shows the dynamic contact stress distributions for horizontal 
vibration and rocking. In all cases, the preceding comment about the 
singularities at the edge holds true. 
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(b) Dynamic Response of a Machine Foundation Embedded 
in the Elastic Halfspace 

In order to show the applicability of the present implementation for 

solving real engineering problems, the dynamic stiffnesses of a machine 

foundation (made of concrete) embedded in the halfspace are computed. 

Dynamic contact stress distributions at the interface between the 

foundation and the soil are also presented. The machine foundation and the 

boundary element discretization for this problem are shown in figure 4.12. 

The discretization of the soil free-surface are the same as in figure 4.3. 

The substructuring technique is used in solving this problem, i.e., the 

concrete foundation is modeled as one GMR (or region) and the halfspace as 

another GMR. The contact between the foundations and the soil is assumed 

to be welded (or glued), and the weight of the foundation is considered 

automatically by the analysis. This problem has corners and edges, and 

therefore, enclosing elements are used to obtain the diagonal blocks of the 

F matrix. The conventional approach of using 0.55„ as the diagonal 

blocks cannot be used for this type of problem which has corners and edges. 

The material properties are as follows: 

Soil: Elastic modulus, E g = 8.64 x 10 5 

Poisson's ratio, \> = o.3 
s 

Mass Density, p = 3.57 

b 

g 

Foundation: Elastic modulus, E c = 4.527 x 10 

Poisson's ratio, v = o.i7 
Mass density, p c = 4.5 

In order to compute the foundation stiffnesses, unit displacements and unit 
rotation are prescribed on the top face of the foundation with zero 
traction conditions being imposed along the soil free-surface. Upon the 
solution of boundary equations, the tractions over the element at the soil- 
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foundation interface are obtained. The resultant of these tractions for 
different rigid bod/’ motions gives the foundation stiffness coefficients. 
The real and imaginary parts of the stiffnesses (minus the inertial 
contribution of the foundation block) are plotted against the frequency 
parameter a Q in figures 4.13 and 4.14, respectively. It can be seen 
that, in general, the stiffnesses in this case are greater than that of a 
rigid strip. This is understandable, because the embedment reduces the 
maximum frequency response (Ref. Estorff and Schmid, 1984) and therefore 
increases the stiffness. Figures 4.15 and 4.16 show the real and imaginary 
parts of the contact stresses between the foundation and the soil for 
vertical vibration whereas figures 4.17 and 4.18 show the same for rocking 
of the foundation. It is obvious from the results, that at higher 
frequency the stresses at the edge are more severe. 

(c) PypamicJ^esROD g e-Of-a Ja il on a n S les tic h alf- spa ce 

Sabjyt gd ^p_a. T im e , Ha nn Qn ig Utera j Pre ggm .e-Pi5tribBtion 

A wall with its base embedded in an elastic half-space is subjected to 
a time-harmonic lateral pressure distribution as depicted in figure 4.19. 
The dimensions of the wall and its base are shown in figure 4.19. The 
material properties of the wall, its base and half-space are the same as 
those of the machine foundation of example (b). 

The boundary element discretization of the wall consists of 20 
quadratic line elements, and its base is modelled by 17 quadratic line 
elements. The discretization of the soil free-surface is the same as in 
figure 4.3. Plane strain conditions are assumed for the present problem. 

The distribution of the applied lateral pressure is shown in figure 
4.19. It can be seen that it is a triangular pressure distribution with 
maximum pressure P(t) =600 psf at the free end of the wall. This 
problem is analyzed by using two different approaches to model the half- 


55 



space, namely, (i) continuum model , and (ii) spring -da shpot model. For 
the spring-dashpot model, the values of stiffness coefficients are 
calculated by assuming the base of the wall to be rigid, and using the 
present dynamic algorithm by following a procedure similar to that 
described in example (a). The lateral displacements along the loaded face 
of the wall are plotted in figure 4.20. From this figure, it can be seen 
that the results obtained by using spring-dashpot model are almost similar 
to those obtained by using continuum model for the half-space. This 
example shows the usefulness of the present algorithm for obtaining the 
response of a structure partially embedded in a half-space in one single 
step or in two-steps, i.e. by using spring-dashpot approach. 

In all of the examples presented in this section, the material damping 
is neglected because for halfspace problems the radiation damping is 
dominant and the material damping is negligible. However, the present 
implementation has the capability for the inclusion of material damping 
(Ref. Sec. IV. 4). 

IV. 7 CONCLUDE REMARKS 

An advanced implementation of the direct boundary element method for 
dynamic analysis of two-dimensional problems in the frequency domain is 
presented. By comparing the results with those obtained by other methods, 
the accuracy and the stability of the present method is established. Since 
only the boundary of the region of interest has to be discretized instead 
of the whole domain, the proposed methodology is a better alternative to 
the conventional finite element method, particularly for the solution of 
soil-structure interaction problems. For soil-structure interaction 
problems the finite element method presents two restraints: (i) the model must 
be bounded at the bottom by rigid bedrock, and (ii) the soil away from 
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the vicinity of the foundation is represented by parallel layers unbounded 
on the horizontal direction. These two conditions are not always close to 
reality. On the other hand, m Boundary element method, the fundamental 
solution satisfies the radiation condition at infinity and therefore no 
bounding surfaces are needed and only a small number of elements is 
necessary to model the problem. Furthermore, the numerical implementation 
employed here is one of the most general currently available and can be 
used in conjunction with substructuring technique to solve not only the 
problems of layered media and soil-structure interactions, but also any 
problem of two-dimensional solids of complicated geometry and connectivity. 


57 



Table 4.1. Vertical compliance of a rigid strip footing on half-space, 

by using three different methods to obtain the diagonal 
blocks of [F] matrix 



Real part of 
the stiffness 
at 

a„ = 2.0 
0 

Imaginary part 
of the stiffness 
at 

a Q - 2.0 

Real part of 
the stiffness 
at 

a Q = 7.0 

Imaginary part 
of the stiffness 
at 

a = 7.0 
0 

Method l 

0.330 

2.24 

0.408 

7.85 

Method 2 

0.335 

2.27 

0.456 

7.81 

Method 3 

0.334 

2.27 

0.410 

7.81 


Method 1: 
Method 2: 
Method 3: 


using 0.5 6^ + j (p_ - F?..)ds as the diagonal blocks. 
AS 

using 0.5 6.. as the diagonal block, 
using enclosing elements. 
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CHAPTER V 

FREE VIBRATION ANALYSIS OF TWO-DIMENSIONAL PROBLEMS 
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V. 1 INTRODUCTION 


In this chapter a new method for free-vibration analysis by BEM is 
presented. It utilizes a fictitious vector function to approximate the 
inertia forces and then uses the well known concept of complementary 
functions and particular integrals to solve the resulting governing 
differential equations. This method not only reduces the problem of free- 
vibration to an algebraic eigenvalue problem but also saves the computation 
time by having fewer matrix manipulations as compared to that of the domain 
integral transformed method (outlined earlier in III. 4). Because of the 
generalized form proposed here it can be used for multi-region problems and 
extensions to axi-symmetr ic problems as well as those involving 
inhomogeneity and anisotropy are possible. Some example problems, such as 
a triangular cantilever plate, a square cantilever plate, a cantilever 
beam, a shear wall and a fixed elliptic arch are presented to establish the 
accuracy, efficiency and convergence of this new method. 


V.2 GOVERNING EQUATION : 

free-vibration of an elastic, 
as: 

( 5 . 1 ) 

where: 

X and 4 are Lame's constants, 
u^ = displacement amplitudes 
p = mass density 
o) = natural circular frequency. 


The governing differential equation for 
homogeneous and isotropic body can be written 


(x + p) 


,2 2 
3 u i 3 “i 2 

5^ * 11 ♦ »“ “i - ° 
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V.3 PARTTCTTTAR INTEGRAL : 


The governing differential equation (5.1) can also be written in 
differential operator notation as 

Lfu^ + pw 2 ^ = 0 (5.2) 

The solution of the above equation can be represented as the sum of a 
complementary function u? satisfying 

L(U?) = 0 (5.3) 

and a particular integral u? satisfying 

L(u?) + p b ) 2 u i = 0 (5.4) 

However, equation (5.4) still contains the unknown displacement field 
within the domain, which can be eliminated by using an unknown 
fictitious density function 6 and a known function C , exactly as in an 
indirect boundary element analysis (Ref. Banerjee and Butterfield, 1981). 
More specifically: 

co 

= 1 C ik (2£.i m )\a m > (5 * 5) 

m=i 

where is a fictitious density and 

C ik is a known function which can be selected as any linear 
function of spatial coordinates. 

The above approximation in the inertia term is a valid practice in other 
numerical methods such as the use of lumped mass matrix in finite element 
method. This is possible because the inertia term does not contain any 
derivative and, hence, it can be approximated by using simpler functions. 
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in the present analysis 


A simple function which is selected for 
has the form: 

c ik<S’i”> - 5 ik« - r) 
where: 

R - largest distance between two points of the body 

r » distance between x (field point) and (source point). 

On the basis of above assumption (5.5), equation (5.4) can be written 
as 


L(u?) + po 2 J C ik (x.I m )^ ]c (i m ) = 0 (5<7) 

m=l 

Now, the particular integral u? can be chosen as any function which 
satisfies the differential equation (5.7). Accordingly it can be 
represented as: 

m 

uP(x) = J (5 - 8) 

m=l 

The displacement field satisfying equation (5.7) is found to be 

pu 2 

D ik (x.i m ) = — [(Cjt - CjRU^r 2 - c^r] 

where: ^i = x i ” ^i 

2(d+3) (l-v)-l 

c = 

18(3d-l) (1-v) 
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JL 


1-2 v 


C 2 = 


2 [ (l+d)-2vd] 


C 3 = 


6(3d-l)(l-v) 


y = Poisson's ratio, and 


d = dimensionality of the problem (e.g. for 2D problems, d = 2). 


For 2-D analysis; 


D ik ■ T [ 'Htp-S ' - 8 ik' 2 - wl <5 - 9 > 


By comparing the functions C^ k and D^ k (for 2D) with the corresponding 
functions f m (eq. 3.6) and (eq. 3.10) of Nardini's method, it can be 
seen that even though the functions and f m are similar, their 
displacement functions D^ k and are different from each other. One 
of the reasons for this difference is that the function D^ satisfies the 
governing differential equation (5.7) but the function does not. 
Instead, the function satisfies the differential equation (3.8) which 
has the form: 

L(^i) = 5 li f Itl (5.10) 

The surface traction t^ related to the displacement u? can be 
determined using the strain-displacement relationship and constitutive 
equation and is given by: 


t?(x) = 2 T ik (x,I m )d k (i m ) (5.11) 

m=l 
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where : 


T ik (x.I m > = po> 2 [<c 4 r - c 5 R)y k n. + (c g r - 2C 2 R) yi n k 
+ ((C 6 c - 2c 2 R)5 ik - 2c 3 y 1 y k /r)y j n j ] 


and: 

(d+3) V - 1 

C 4 = 

4 3(3d-l)(l-v) 


2 v 

(i+d)-2vd 

(d+2)-(d+3)v 
3(3d-l) (1-v) 


(5.12) 


V.4 R Of IMPART FT.EMEyrT FORMULATION : 

The boundary values of real displacements and tractions u^ and t^ 
can be related to the complementary and particular integral via: 


“i-"'* 4 

(5.13) 

tl - t= + tf 

(5.14) 


The boundary integral equation related to the displacement function 
u? can be written as 

c,.(I)u?(l) « f [g. .( x.i)t?(5) - F, .(x.I)u?(x)] ds(x) (5.15) 

Ij 1 J U lj * ij * J 

s 

where G-^tx.jL) and are the fundamental solution of equation 
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(5.3) (Ref. Barter jee and Butterfield, 1981), i.e., G^(x,£) is the 

resulting displacement at any point £ in direction i of an infinite 
medium due to a static unit force acting at a point I in direction j , 
and F^ (x,i) is the resulting traction. 

By usual discretization of boundary S , we can express equation 
(5.15) in matrix formas 

(GHt c ) - [FHu c } = CO) (5.16) 

Equation (5.16) can be solved once the displacements u*r and the 
tractions t^ are expressed using equation (5.13) and (5.14) in terms of 
real displacement u^ and traction t^ » i.e. 

CGI Ct) - [FHuJ = [GHtP) - [FHu p ) (5.17) 

where vectors Ct p } and tu p } can be obtained at boundary nodes from 
equations (5.8) and (5.11) as 

(u p ) = pw 2 [DHdl (5.18) 

{t p } = p« 2 CT] Id) (5.19) 

Substituting these equations into equation (5.17), we obtain 

[Gift) - tFHu) = pu 2 ( [G] [T] - CFl CD] ) Cd) (5.20) 

Recalling that 

<D 

u i(x n ) = J 6 t j (R - r 1 ™)^^) 
m=l 

where r™ 11 is the distance between the points * n and £. m , we can 
express this relationship between the displacements and the fictitious 
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density at all boundary nodes as: 


u i = S ij pnm<z5 ' j 1 

d? = 8 ij Rnmu i (5.21) 

where K m = (P™) -1 

It should be mentioned here that we only have to invert a NxN matrix 
(Pi instead of [Ql (as in the case of Nardini's method, eq. 3.12), a 
dNxdN matrix, where d = 2 and 3, respectively for two and three- 
dimensional problems and N is the total number of boundary nodes of the 
problem. 

We can now write (5.21) as 

M) « [KHu) (5.22) 

Substituting (d) from equation (5.22) into equation (5.20), we get 

[GHtJ - [FHu] = pa» 2 ([G][T] - [F] [D] ) [K] tu) (5.23) 

or IGHt} - [FHu] = p<o 2 [M] (u) (5.24) 

Equation (5.24) can also be written in terms of known and unknown variables 
as 

lAHx) - IB] (y) = pu 2 ( [M] (x) - [M*Hy}) (5.25) 

Since all the known variables are zero, (i.e. specified boundary conditions 
are either u^ = 0 or t^ = 0 ) equation (5.25) reduced to 


The modified mass matrix [M] contains zero in its sub-columns related 


to specified displacements (i.e. fixed boundaries). 


V. 5 EIGENVALUE EXTRACTION : 

Equation (5.26) is an algebraic expression for the eigenvalue problem 
which can be solved by using a eigenvalue extraction subroutine. It should 
be noted that both the matrices [A] and [M] are fully populated and 
nonsymmetric. There is no satisfactory eigenvalue extraction routine 
available for efficient determination of eigenvalues of such a system. In 
the present work the algorithm developed by Moller and Stewart (1973) was 
utilized. The necessary set of subroutines were developed by Garbow (1980) 
of Argonne National Laboratory. In general a nonsymmetric fully populated 
system such as (5.26) cannot be guaranteed to provide real eigenvalues. 
However, it will be seen from the examples presented in this chapter that 
the eigenvalues of (5.26) are in fact real. 

V.6 ADVANTAGES OF THE PROPOSED METHOD : 

In comparing this new method with that of Nardini's (1982), the 
following three important points need to be mentioned: 

(i) The final algebraic expression of Nardini’s method is in terms of 
unknown displacements whereas that of this new method is in terms of 
unknown variables (both displacements and tractions). In a multi- 
region (piecewise homogeneous) problem both the displacements and 
tractions are unknown at the interface. Therefore, Nardini's 
assumption that at ary node either the displacement or the traction is 
zero is not always valid. 

(ii) Nardini's approach involves too many matrix manipulations which are 
costly and somewhat impractical for a realistic practical problem. 
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(iii) Because of the use of a well-established method of solving any 
nonhomogeneous differential equation by using complementary functions 
and particular integrals, we can now utilize a large number of 
particular integrals already in use in BEM for dealing with 
centrifugal forces. Many of these have already been developed for 
axi-syrrmetric and three-dimensional problems involving anisotropic 
materials and, with minor modifications, can be made to satisfy the 
governing differential equation (5.7). 

V.7 EXAMPLES OF APPLICATIONS : 

(a) Comparison with Nardini and Brebbia (1982) 

In order to provide a meaningful comparison between the present method 
and that described by Nardini and Brebbia, both methods were implemented. 
Tables 5.1 and 5.2 show the convergence studies of the first four modes of 
triangular and square cantilever plates of unit thickness under in plane 
vibration. The triangular plate had a 10-inch depth at the support and an 
8-inch span. The square plate was 6 inch deep and had a span of 6 inches. 
The material parameters were E/p= 10 4 and \> = 0.2. Three-noded, 
isoparametric-conforming boundary elements were used to describe both 
geometry and functions. 

Both these problems were also solved by Nardini and Brebbia (1982). 
The results of the present implementation agree exactly with their quoted 
results indicating that their analysis has been correctly interpreted. 
They do not, however, agree well with those given by the new method 
proposed in this paper for some modes. Specifically, for the triangular 
cantilever, there is a marked difference in the third mode and small 
differences exist in all other modes. For the rectangular cantilever once 
again third-mode response differs significantly but the second mode agrees 


68 


quite well. 


(b) Comparisons with Finite Element and Beam Theory 

The finite element system MHOST (MARC-HOST) was used to analyze a 
cantilever beam. The beam has a length of 6.5 units and a square (1 x 1) 
cross section. The finite element mesh (using 8-noded isoparametric 
elements) was matched with the boundary element mesh to provide the same 
number of boundary nodes. The first four bending modes from BEM were 
(0.368. 2.214. 5.591 and 9.986 Hz) and those of the FEM were (0.378. 2.188, 
5.583 and 9.908 Hz), indicating good agreement between the two analyses. 

Further the mode shapes calculated using the two techniques are 
indistinguishable. The first and the fourth bending modes are shown in 
Figure 5.1. It should be noted that the fourth mode displays a nonzero 
slope near the fixed end. This real feature of the two-dimensional 
solution is absent in the beam theory with the imposed fixed end boundary 
conditions normally used in the beam theory. The material parameters for 
the beam are assumed to be E/p = 10 4 and v = 0 . 

In order to study the convergence of the results with an increase in 
number of boundary elements, a similar cantilever having a span of 6.0 
inches was analyzed. Figure 5.2 shows the convergence of the first six 
modes plotted against the boundary mesh numbers. Total number of boundary 
elements is equal to 2x(Mesh number + 2) (Ref. Fig. 5.3). The convergence 
is excellent for the first six modes. Since this analysis is fully two- 
dimensional rather than based on beam or column theory, it provides both 
the axial and flexural modes. In addition, seme of the higher modes (not 
shown here) have mixed responses. As expected, a finer discretization is 
required for higher modes of vibration. Even the most slowly convergent 
case, the fifth (the fourth flexural) mode required only 8 boundary 
elements. This indicates that the present analysis could be further 
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developed to provide a powerful analytical tool for free vibration 
analysis. 

Figure 5.2 also shows the beam theory results for modes 1. 2 and 4 
(the first three flexural modes). The increasing departure of the results 
from the beam theory is due almost entirely to the neglect in beam theory 
of shear distortion. Approximate modifications of the beam theory results 
for a simply supported beam (Ref. Clough and Penzien, 1975) to account for 
this effect indicate frequency reductions of approximately the magnitude 
observed. 

(c) An Example of a Shear Wall 

In order to compare the results obtained from the proposed method with 
those from the Finite element method and Nardini's BEM, a shear wall with 
four square openings was analyzed for in-plane vibration. The boundary 
element and the finite element meshes (Ref. Nardini and Brebbia, 1982) are 
shown in Figure 5.4. The material parameters were E/p = 10 4 and 
v = 0.2. 

Table 5.3 shows free- vibration periods for the first eight modes. The 
first mode is identically same as that obtained by FEM. The present 
results for 2nd. 3rd and 5th modes are also close to the FEM results. The 
results from the present analysis agree well with those reported by Nardini 
for 4th. 6th, 7th and 8th modes. However, they do not agree well for the 
rest of the modes. 

(d) An Example of an Arch with Square Openings 

An arch in plane stress and fully fixed at the supports was analyzed 
(Figure 5.5) for in-plane vibration. Four different cases involving the 
full arch with or without openings and symmetric halves with or without 

n 

openings were considered. The material parameters were E/p = 10' and 
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V = 0.2 . 

Table 5.4 shows the natural frequencies of the full arch with and 
without openings. In general the natural frequencies are reduced due to 
the presence of openings which affects some modes more than others. 
Similar results for the symmetrical half of the arch are shown in Table 
5.5. In this latter case, of course, some of the nonsymmetric modes of the 
full arch are absent. Modes 1 and 6 of the full arch are identical to the 
first two modes of the symmetric half. 

V.8 CONCLUDING REMARKS : 

A new method based on the well known technique of solving a 
nonhomogeneous differential equation by complementary function and 
particular integrals for the analysis of free vibration problems by 
boundary element is presented. The method has been compared with that of 
Nardini and Brebbia (1982) and found to yield different results for some of 
the higher modes of vibration. It has also been compared with MARC-HOST 
finite element analysis and was found to yield essentially similar results 
for a cantilever beam problem. When the beam theory is corrected for the 
shear deformation, the analytical results tend to agree well with those of 
the present analysis. 

The present analysis can be easily extended to axi-symmetric and 
three-dimensional problems involving inhomogeneity and anisotropy by 
utilizing a number of particular integrals already in use in boundary 
element analysis. 
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TABLE 5.1: Time Periods of 

Number Mode 1 


of 

Slements 

Nardini's 

method 

New 

method 

3 

0.415 

m 

6 

0.415 


9 

0.416 

0.430 

12 

0.416 

0.430 

15 

0.416 

0.430 

18 

— 

0.430 



TABLE 5.2: Time Periods of Free Vibration of a Square Cantilever Plate 


Number 

Mode 1 

Mode 2 

Mode 3 

Mode 4 

of 

Nardini's 

New 

Nardini's 

New 

Nardini's 

New 

Nardini' 

s New 

Elements 

method 

method 

method 

method 

method 

method 

method 

method 

4 

0.536 

0.561 

0.232 

0.235 

0.195 

0.172 

0.109 

0.107 

6 

0.545 

0.568 

0.234 

0.237 

0.214 

0.179 

0.118 

0.116 

8 

0.559 

0.581 

0.236 

0.238 

0.210 

0.185 

0.127 

0.122 

10 

0.562 

0.581 

0.236 

0.238 

0.209 

0.187 

0.129 

0.123 

12 

0.563 

0.584 

0.237 

0.238 

0.209 

0.187 

0.131 

0.125 

16 

— 

0.585 

— 

0.238 

— 

0.187 

— 

0.125 


TABLE 5.3: Time Periods of Free Vibration of a Shear Wall 


Modes 

1 

2 

3 

4 

5 

6 

7 

8 

FEM 

(SAPIV) 

3.029 

0.885 

0.824 

0.526 

0.409 

0.342 

0.316 

0.283 

Nardini's 

BEM 

3.022 

0.875 

0.822 

0.531 

0.394 

0.337 

0.310 

0.276 

New 

Method 

3.029 

0.878 

0.823 

0.533 

0.400 

0.337 

0.311 

0.276 



TABLE 5.4: Free vibration modes of full arch without and 

with openings (Hz) 


Modes 

Without openings 

With openings 

1 

87.8 

78.9 

2 

124.1 

113.5 

3 

177.4 

146.8 

4 

230.9 

212.4 

5 

275.7 

235.0 

6 

380.7 

265.5 

7 

428.1 

401.1 

8 

506.1 

537.1 

9 

622.0 

590.9 

10 

648.0 

595.2 


TABLE 5.5: 

Free vibration modes of the symnetric half of the arch 
without and with openings (Hz) 

Modes 

Without openings 

With openings 

1 

123.9 

113.4 

2 

378.9 

264.4 

3 

429.3 

395.7 

4 

649.4 

590.5 

5 

820.1 

670.6 
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CHAPTER VI 

ADVANCED THREE-DIMENSIONAL STEADY-STATE DYNAMIC ANALYSIS 
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VI. 1 INTRODUCTION 


In this chapter, an advanced implementation of the direct boundary 
element method applicable to the steady-state dynamic analysis of problems 
involving three-dimensional solids of arbitrary shape and connectivity is 
presented. Isoparametric curvilinear surface elements are used for mapping 
geometry and for approximating variation of the field variables. In the 
present implementation, substructuring capability is incorporated for 
solving problems involving piecewise-homogeneous materials such as problems 
of layered media and soil-structure interaction. Also provided is a 
feature called built-in-symmetry; this allows one to solve the problems 
having geometric and loading symmetry by modelling only a part of the 
actual geometry. In this chapter, the discussion starts with the boundary 
element formulation for steady-state dynamics followed by techniques 
related to the numerical implementation. The assembly and solution 
algorithms for general three-dimensional problems are the same as those for 
two-dimensional problems (Ref. Secs. IV.4.G and IV.4.H), and therefore they 
are not repeated in this chapter. Finally, a number of numerical examples 
are presented to demonstrate the accuracy and applicability of the present 
implementation. This dynamic analysis technique seems to provide an 
accurate and efficient tool for solving truly three-dimensional problems 
and particularly those relevant for problems of soil-structure interaction, 
where it has clear advantages over existing finite element solutions. 

VI. 2 BOUNDARY INTEGRAL FORMULATION 

The boundary integral equation for three-dimensional problems of 
steady-state elastodynamics is the same as that of two-dimensional problems 
(eq. 4.13) and it can be expressed as: 
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( 6 . 1 ) 


Cij (i)Ui<i,M) = J CG^j (x»I.“)u i (S.“)]<3S(x) 

s 

The above equation is identical to equation (4.13), except that here 
the field variables and the fundamental solution are functions of circular 
frequency w rather than that of the Laplace parameter s . This is 
possible because s and <■> are interchangeable (s = -iw) . The 
fundamental solution G^j and F — are listed in Appendix A2. It should be 
noted here that although the functions G^ and F — becomes identical to 
their static counterpart as s tends to zero, it is important to evaluate 
this limit carefully because of the presence of s in the denominator. 

Once the boundary solution is obtained, the stresses at the boundary 
nodes can be calculated by combining the constitutive equations, the 
directional derivatives of the displacement vector and the values of the 
field variables at the boundary nodes in an accurate matrix formulation 
(Ref. Sec. VI. 3. G). Also the loads and moments can be obtained by 

numerically integrating the known tractions on each element. 

For displacements at interior points, equation ( 6 . 1 ) can be used with 
appropriate tensor (Ref. Sec. IV.3); and the interior stresses can be 

obtained from 

= j [G?j k (x,JL,w)t^(x.«>) - F^j k (x»i.“)u^(x,(o)]dS(x) 

( 6 . 2 ) 

The functions G?^ and F^^ in the above equation are listed in 
Appendix A3. 

The constitutive equations and boundary conditions are the same as 
described in Chapter IV (Ref. Secs. IV.1-IV.3). This boundary integral 
formulation presented above can also take account of viscous damping (Ref. 
IV.3). 
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IV. 3 NUMERICAL IMPLEMENTATION ; 

Since the basic governing equations for dynamic analysis in the 
transformed space (either in w or s space) are similar to the 
corresponding equations for the static analysis, the numerical 
implementation developed for the static case can be used to extract 
solution for the dynamic problem for one value of the Laplace transform 
parameter s or frequency parameter w . The current numerical 
implementation of the boundary integral equation for three-dimensional 
problems of steady-state dynamics has the following aspects and features. 

(A) Representation of Geometry and Field Variables 

The boundary integral equation (6.1) represents an exact formulation 
involving integrations over the surface of the domain. Therefore, if one 
does not make grossly simplified assumptions in the spatial variations of 
the boundary quantities, accurate solutions can be obtained. To this 
purpose, each surface is discretized in a number of elements with each 
element defined in terms of several geometric nodes. All surface-element 
types employed represent surface geometry using quadratic shape functions. 
Three sided elements, defined using six rather than eight geometric nodes, 
are used for mesh transition purposes. The terms quadrilateral and 
triangle are normally used to refer to the eight and six noded elements 
although the real geometry represented is, in general, a nonplanar surface 
patch in three dimensions (Ref. Fig. 6.1). Over each element the variation 
of field variables can be defined using either the linear or quadratic 
shape functions. Linear and quadratic elements can share a common side 
which is then constrained to have linear displacement and traction 
variation. 

In addition to the element types mentioned above, elements which 
extend to infinity are provided. These elements are designed to allow 
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modeling of structures connected to the ground and automatically 
incorporate appropriate decay conditions. The characteristics of the 
various element types are summarized below (Ref. Appendix C2). 


Element type G eom et ry .; N odes Field variable Nodes 


Linear Quadrilateral 8 4 
Linear Triangle 6 3 
Quadratic Quadrilateral 8 8 
Quadratic Triangle 6 6 
Quadratic Infinite 8 3 


The cartesian coordinates x-^ of an arbitrary point P on a surface 

element are given in terms of the nodal coordinates X. as: 

la 

x i ( p) = N a O.)X ia (6.3) 

where i = 1,2,3 and a = 1,2,. ...A, with A the number of nodal points 
necessary to describe the element. Furthermore, N a are the shape 
functions defined in the local or intrinsic coordinate system The 

Jacobian matrix relating the transformation from the cartesian coordinate 
system (x,y,z) to the element's intrinsic coordinate system (t^.t^) is 



(wr 0 /*n j )x itt 


(6.4) 


where j - 1,2 and the summation convention is again implied for repeated 
indices such as a . 

The field variables are also represented by the same shape functions, 

i.e. 


u. (x) = N (ii)U- and 

1 — a -*• la 

t t (x) = N a (a)T ia (6.5) 
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where U^ a and T^ a are nodal values of the displacements and tractions, 
respectively, in the transformed domain. 

Infinite elements, which are essential if problems involving the half- 
space are to be solved, can be constructed by modifying the eight-node 
quadrilateral as shown in Fig. 6.3. The intrinsic coordinate along the 
dimension of the quadrilateral that we want extended to infinity (say ri 1 ) 
is modified as 

nj = (3n x + i)/(i - u 1 ) (6.6) 

This way the original interval (-1,1) is mapped into (-1,+®). It should be 
noted that only the three nodes on the side of the infinite element that is 
adjacent to a surface element belonging to the 'core' region contribute to 
the system equations. The original shape functions N for these three 

nodes are then modified by the ratio d = [(x^-z ^)(xj-z^)/(y^-z^) (y^ - 

1/2 2 
z^)] for the displacement kernel and d for the traction kernel where 

x^ are the cartesian coordinates of the integration points, y^ their 

projection on the common side with the core, and z^ an arbitrary 

reference point. This type of stretching of a quadrilateral results in a 

Jacobian determinant equal to 4/(l-n 1 ) 2 that must be included in the 

kernel integrations. The infinite element thus obtained reproduces the 

correct spatial decay of the fundamental singular solutions as r -> » . 

(B) Built-in Symmetry and Sub-structuring Capabilities 

In obtaining the numerical solutions, the built-in symmetry capability 
allows one to solve the problems having geometric and loading symmetry by 
modeling only a part of the actual geometry. The major steps in this 
procedure are briefly explained as follows. If the geometry and the 
boundary condition are symmetric with respect to a plane (or a number of 
planes), then only that portion of the boundary which lies on the one side 
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of the plane (or planes) is modeled. The symmetry can be with respect to 
y-z plane (half -symmetry), y-z and x-z planes (quadrantal symmetry), or y- 
z, x-z and x-y planes (octan symmetry). The effect of the unmodeled part 
of the boundary is included according to the following scheme: For all the 
field points, the contribution of the unmodeled portion to the matrices of 
coefficients F— and are accounted for by reflecting the modeled 

surface elements with respect to the plane (or planes) of symmetry and then 
integrating over the reflected elements (with proper normals). For the 
source nodes on the plane (or planes) of symmetry the contributions are 
added up directly whereas for all other source nodes the correct signs of 
the contributions are determined by the directions associated with the 
field variables with respect to the plane (or planes) of symmetry. By 
avoiding the calculation of identical quantities, this procedure shortens 
the time required to evaluate the matrices. In addition, it reduces the 
time required to solve the set of linear equations, because the system 
matrix will have fewer rows and columns. 

The substructuring capability allows a structure to be modeled as an 
assembly of several generic modeling regions (GMR). The GMRs, each of 
which must be a complete portion of the structure, are joined by enforcing 
appropriate compatibility conditions across common surface patches 
(elements). This feature can also be used to solve piecewise inhomogeneous 
problems because the GMRs can have different material properties. 

(C) frfc mecic al In tegra t io n 

In view of the surface elements introduced in the previous section, 

Eq. (6.1), when integrated over the surface of the problem in question, assumes 
the following form: 
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Q 

c ij (i , “ i ( i) = } { I 6 ij (x(a).I.s)N 0 (a)ds(x(n)TT o 

< 3= 1 S q 

- f F^(x(a).i.s)N (aJdSCxfa))^. } (6.7) 

J ij u lu 

S q 

In the above equation. S is the surface of the qth element and Q is the 
total number of elements. The global system of boundary element equations 
at a given value of s is obtained fcy the usual nodal collocation scheme, 
i.e.. by allowing point i to coincide sequentially with all the nodal 
points of the boundary. 

With the exception of strongly singular traction integrals, all 
surface integrals in the numerical implementation have been calculated 
numerically. Since this is the most time consuming portion of the 
analysis, it is essential to optimize this effort. Essentially two types 
of integrals, singular and nonsingular, are involved. Hie integrals are 
singular if the field point for the equations being constructed lies on the 
element being integrated. Otherwise, the integrals are nonsingular 
although numerical evaluation is still difficult if the field point and the 
element being integrated are close together. 

In both the singular and nonsingular cases, Gaussian integration is 
used. The basic technique is developed in Banerjee and Butterfield (1981) 
and was first applied in the three-dimensional boundary el orient method by 
Lachat and Watson (1976). in the nonsingular case an approximate error 
estimate for the integrals was developed based on the work of Stroud and 
Secrest (1966). This allows the determination of element subdivisions and 
orders of Gaussian integration which will retain a consistent level of 
error throughout the structure. Numerical tests have shown that the use of 
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3x3 , 4x4, and 5x5 Gaussian rules provide the best combination of accuracy 
and efficiency. In the present implementation the 4x4 rule is used for 
nonsingular integration and error is controlled through element 
subdivision. Typical element subdivisions into three-node triangles and 
four node quadrilaterals are shown in Fig. 6.4. The distance R that 
controls the subdivision process is measured frcm the field point to the 
point closest to the field point on the element being integrated. In 
general, higher values of s require lower integration tolerance leading 
to more element subdivision. If the field point is very close to the 
element being integrated, use of a uniform subdivision of the element can 
lead to excessive computing time. In order to improve efficiency while 
still retaining accuracy, a graded element subdivision is employed. Based 
on one-dimensional tests, it was found that the subelement divisions could 
be allowed to grow geometrically away from the origin of the element 
subdivision. Numerical tests on a complex three-dimensional problem have 
shown that a mesh expansion factor as high as 4.0 can be employed without 
significant degradation of accuracy. 

In the case of singular integration, which arises when the field point 
is on the element being integrated, the element is first divided into 
triangular sub-elements. The integration over each sub-element is carried 
out in a polar coordinate system with the origin at the field point. This 
coordinate transformation produces nonsingular behavior in all except one 
of the required integrals. Normal Gaussian rules can then be employed. 
The integral of the traction kernel times the isoparametric shape function 
which is l.o at the source point is still singular and cannot be 
numerically evaluated with reasonable efficiency and accuracy. Its 
calculation is carried out indirectly as discussed in Chapter IV, Section 
4JD. It has been found that subdivision in the circumferential (angular) 
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direction is required to preserve accuracy in the singular integration. A 
maximum included angle of 15 degrees is used. Subdivision in the radial 
direction has not been found necessary. This process is illustrated in 
Fig. 6.5 for a quadrilateral element. 

The surface integrals required for calculation of displacement and 
stress at interior points are of the same type as those involved in the 
boundary problem with the exception that only nonsingular integrals are 
involved. In general, the integrals appearing in the surface integrals are 
continuously differentiable and solution accuracy can, therefore, be 
improved by use of increased integration order. 

(D) Calculation of Stresses on the Boundary for 3D Problems 

Once the boundary solution is obtained, the stress and strain at any 
point on the boundary can be calculated without any integration, by using 
the procedure outlined as follows. 

Let us assume that we are interested in finding stress and strain at a 
point P , which lies in a boundary element and has intrinsic coordinates 
Recalling equations (6.5), we can write 

A 

u i (l »l' T »2 ) = l V t, l* T, 2 )u i« 

a=l 

A 

Vn$.1j> ■ } Vl-’ l 2 )t l« (6 - 8> 

a=l 

where: 

A is the number of nodes in the element, 

N q is the shape functions, and 

u. and t • „ are the nodal values of u< and t^ . 
la la 11 

In addition, we also have the following relationships: 


85 



(6.9) 


fc i = CT ij n j 


_ U k,l +U l,k 
"ij - C ijkl ( 2 } 


( 6 . 10 ) 


u i.n =u i.j X j.n 


( 6 . 11 ) 


u i., ■ l 


3N 
a 

3tj U io 


( 6 . 12 ) 


where 

c i j )cl is a tensor containing elastic constants, and 
X. „ are the directional derivatives. 

Equations (6.9), (6.10) and (6.11) can be combined to form a matrix 

equation: 


[Slip) = (q) 


(6.13) 


where IS] is a 15x15 matrix which contains unit normals, a 3x3 unit 
matrix and material constants; {p} is the unknown vector of <r^ and 
3u^/ 3 5 j ; and (q) is a vector containing the tractions t^ and local 
derivatives of the displacements at point P . 

Finally, the stress and strain at point P can be obtained by 
inverting the matrix of equation (6.13) and then multiplying the inverted 
matrix by the right-hand-side vector. For this purpose, the right hand 
side vector is obtained by using equations (6.8) and (6.12). 


A number of representative problems were solved in order to test the 
steady-state solution. In all cases, English units are used with foot (ft) 
for length, pound (lbf) for mass, and second (s) for time. 
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(a) Cantilever Subjected to End Shear 


A uniform team with a rectangular cross-section is completely fixed at 

iOT 

one end and a uniformly distributed traction t^ = 1000 e . with n = 314 
r/s, is applied at the other end. Traction-free conditions hold along the 
sides. The dimensions of the beam are L - 10 , w = 1 , and d = 3 . The 

n 

material properties are as follows: modulus of elasticity E = 1.16 x 10 

and mass density p = 2.0 . in order to reproduce the one-dimensional 
characteristics the Poisson's ratio is assumed to be equal to zero. This 
cantilevered beam is modelled by 18 quadrilateral surface elements 
resulting in 56 nodes. In reference to Fig. 6.6, it is observed that the 
surface elements are arranged closer to the loaded end. This is so because 
the displacement function varies more sharply at the loaded end than at the 
fixed end. The same figure plots the absolute value of the vertical 
displacement u y along the length of the beam at a frequency <o equal to 
the forcing frequency n . The results are in very good agreement with the 
analytical solution for a flexural beam which was developed from Clough and 
Penzien (1975). 

(b) Cantilever Subjected to Harmonic Transverse Load 

The same model discussed in (a) was subjected to a time harmonic patch 
load as shown in Fig. 6.7. The agreement between the three-dimensional 
calculation and beam theory (Clough and Penzien, 1975) was, once again, 
excellent. 

(c) vertical Compliance of_a_ Rigid Square. Footing 

A rigid square foundation of side length 2b = 2 is resting on the 
surface of a homogeneous halfspace under relaxed boundary conditions (i.e., 
there is no friction between soil and foundation). The halfspace has a 
shear modulus p = 1.0 , v=l/3, and p = 1.0 . The foundation is 
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subjected to a uniform harmonic vertical displacement u of amplitude 
equal to unity. The surface of the halfspace is traction free. The 
traction distribution under the foundation obtained by the BEM is 
integrated to give the total vertical load P_ • The foundation's 
normalized compliance in the vertical direction is obtained as C vv = 
fibu / P . The two meshes shown in Fig. 6.8 are used for modelling the 
foundation as well as surface of the halfspace. Since the transformed 
domain BEM computer program can take advantage of symmetry, only 1/4 of the 
problem needs to be discretized. The coarse mesh uses 4 and 12 elements to 
model the foundation and the halfspace, respectively. Note that the 
outermost 4 elements are infinite elements. This discretization results in 
44 nodes. The finer mesh uses 6 and 12 elements for the same purpose. 
There are 2 infinite elenents here and 65 nodes. 

This problem was originally solved by Wong and Luco (1976). They 
numerically integrated the vertical displacement at the surface of a 
homogeneous halfspace due to a unit point load over the foundation, which 
was discretized into small squares. This problem was recently revisited by 
Rizzo et al (1985) using a BEM approach. In their work (Rizzo et al), both 
frictionless and welded cases are considered and two approaches are used: 
The exact one employs the halfspace kernels (Lamb's solution) and the 
approximate one uses the fullspace kernels (Stoke's solution). In both 
cases, only the rigid foundation is discretized and these two approaches 
are practically indistinguishable except at the very low frequency range. 
All three solutions mentioned are plotted in Fig. 6.9, along with the 
vertical compliance obtained by the present method using the fine mesh. 
The good agreement between the present results and that of Rizzo et al 
(1985) should be noticed. However, the major difference between Wong and 
Luco's results and the boundary element results is due to the fact that 
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quadratic shape functions are used for representation of the variation in 
the field variables over each element in the present work as well as that 
of Rizzo et al whereas Wong and Luco assumed that the unknown contact 
stresses are uniform within each element which is a crude approximation. 
Finally, the difference in results obtained by both coarse and fine meshes 
is contrasted in Table 6.1. 

VI. 5 CONCLUDING REMARKS 

An advanced algorithm based on the direct boundary element method for 
the steady-state dynamic analysis of structures behaving elastically or 
viscoelastically has been presented. The numerical implementation employed 
is one of the most general presently available and can be used in 
conjunction with substructuring to treat three-dimensional solids of 
complicated geometry and connectivity. The algorithm is stable and capable 
of producing very accurate results except perhaps at high frequencies in 
which case finer meshes are required for better accuracy. Nevertheless, 
the present method is a viable alternative to algorithms based on finite 
element methodology. Specifically for halfspace problems, the present 
method does not require discretization of the domain of the halfspace and 
the use of energy absorbing elements as is required by the finite element 
method. 

The present method can very easily be extended to solve time-harmonic 
wave scattering problems by simply adding the displacements due to the 
incident field on the right hand side of the final system equation. 
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Table 6.1: 

Comparison of vertical compliances obtained by using 
two different meshes 

wb 

Re{c W (a o } 


x®«W 

a o )} 

o' =2 ' 

Coarse Mesh 

Fine Mesh 

Coarse Mesh 

Fine Mesh 

0.5 

0.118 

0.117 

-0.057 

-0.058 

1.0 

0.064 

0.069 

-0.083 

-0.081 

1.5 

0.032 

0.034 

-0.076 

-0.070 

2.0 

0.021 

0.018 

-0.059 

-0.052 

2.5 

0.015 

0.015 

-0.052 

-0.048 

3.0 

0.010 

0.012 

-0.036 

-0.037 

3.5 

0.005 

0.006 

-0.035 

-0.032 

4.0 

0.004 

0.004 

-0.027 

-0.027 
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C HR g E E R V TI 

TRANSIENT DYNAMIC ANALYSIS BY LAPLACE TRANSFORM 
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VII. 1 


niaa>ii!i 


UCTION 


In this chapter, an advanced implementation of the transformed domain 
boundary element formulation applicable to transient dynamic problems 
involving two and three-dimensional solids of arbitrary shape and 
connectivity is presented. Using the correspondence principle (Lee, 1955), 
the transient dynamic problem is first solved in the Laplace transform 
space and then time domain solutions are obtained by numerical transform 
inversion. The transformed governing equations and the transformed 
boundary element formulation are presented in Chapter IV (Sec. 4). The 
materials pertaining to the fundamental singular solutions and the 
numerical impl mentation of the boundary integral equation for one value of 
Laplace transform parameter are discussed in Chapter IV (Secs. 4-5) and 
Chapter VI (Secs. 2-4) for two and three-dimensional problems, 
respectively. This chapter starts with a discussion on the Laplace 
transformed equations of elastodynamics followed by numerical inversion of 
Laplace transform. Numerical examples are finally presented and, through 
comparisons with available analytical and numerical results, the stability 
and high accuracy of this dynamic analysis technique are established. 

VII. 2 LAPLACE TRANSFORMED EQUATIONS OF ELA5?TCDYNAMICS 

The governing differential equation of linear elastodynamics in 
Laplace transform domain can be written as: 

(c 2 - c„ 2 )u- • • + c„ 2 u- • • + b- - s 2 u- + U- + sU. = 0 (7.1) 

1 2 /u i,ij 2 ] , ll 2 3 3° jo 

With the assumption of zero initial condition and absence of body 
force, the above equation reduces to: 
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(7.2) 


(C, 


- C 2 2 )U, 


1 , 1 ] 


+ c. u. 


2 j ,ii 


2 ‘ 

- s u . 


= 0 


Since the boundary condition and the constitutive equations do not 
involve time derivatives, their Laplace transforms are simply: 


u i = q^X.s) 

ti = = p i (x-s) 

ff. ■ = p [(c 2 - 2c 2 )u 5.. + C,*(U. • + U • .)] 

13 l 2 m,m lj 2 1,3 j.i 


(7.3) 


Finally, the boundary integral equation in Laplace transform domain 
has the form 


Cij(£) u i(£,s) = J [G i j(x,i,s)t i (x,s) - F^j(x»l.s)u^(x.s)] dS(x> 

S ' (7.4) 

The main advantage of casting the equations in the Laplace transform 
domain is that the equations of motion become elliptic partial differential 
equations, and as such are more amenable to numerical solutions than their 
hyperbolic counterparts in the time domain. The numerical solution of the 
transient elastodynamic problem in the Laplace-transform domain essentially 
consists of a series of solutions to a static-like problem for a number of 
discrete values of the transformed parameter s . The final solution is, 
of course, then obtained by a numerical inversion of the transformed domain 
solutions to the time domain. 


VII. 3 DIRECT LAPLACE TRANSFORM OF BOUNDARY CONDITIONS 

In order to solve equation (7.4), the boundary conditions have to be 
transformed to the Laplace domain. As the input boundary conditions are 
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piecewise linear in time, a numerical scheme is used to transform the 
boundary conditions from time domain to Laplace domain. The formula used 
for this purpose is exact for the forcing functions (i.e. boundary 
conditions in our case) which are piecewise linear in time, and is given by 


N-l 


f(x.a) = ) {AF(e n - e 


n=l 


s AT 


-sT_ -sT -sT -sT i 

n X ) + sAT(P n e n - F n+1 e n 1 


)} 
(7.5) 


where: 


F_ = f (x»T_) = value of f at time T , and 
n n n 


“■ - P n*l - F n 


The above formula is tested for a number of trial functions (such as 
coswt. e , logT, etc.) for N = 20 and N = 50 . The average error for 
N = 50 is 0.5 percent and that for N = 20 is 1.2 percent. Therefore, 
(7.5) can also be used for taking Laplace transform of any arbitrary 
loading function. 


VII. 4 NUMERICAL INVERSION OF TRANSFORM DOMAIN SOLUTION 

After numerically integrating equation (7.4) over the surface and 
imposing known boundary conditions, the final system equations can be 
assanbled to the form 


[AHX) = [BHY) (7.6) 

All expressions in the above equation are dependent on the transform 
parameter s . Therefore, for a transient dynamic problem, the above 
equation is formed and solved for (X) for a spectrum of values of the 
transform parameter. 
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Finally, all that remains to be done is to invert the solutions back 
to the real time domain. In general, transformation from the Laplace 
transform domain back to the time domain by analytical methods is 
impossible except for simple functions. Therefore, numerical evaluation of 
the inverse Laplace transform is imperative. The inverse Laplace transform 
can be defined as 


f(X,T) 


i r Y+i “ - 

2nT J . f( *' S)e 

y-l® 


ST 


ds 


i = 7-1 


(7.7) 


where y ( > 0) is arbitrary but greater than the real part of all the 
singularities of f(&,s) and s is a complex number with Re(s) 1 y > 0. 

The various methods available for numerical inverse Laplace 
transformation may be grouped (Ref. Narayanan, 1982) as follows: (a) 
Interpolation-collocation methods, (b) methods based on expansion of 
orthogonal functions, and (c) methods based on numerical Fourier 
transforms . 

In this work, Durbin's (1974) method is used because of its high 
accuracy (Ref. Hanoi is et al, 1981; and Ahmad et al, 1985). Durbin's 
method is classified under group (c) and combines both the Fourier sine and 
the cosine transforms to arrive at the inversion formula: 


f(x,T.) = 2(^r-)f- 0.5 ReCf(x.y)} + Re { J (A(n) + iB(n))W^ n }] 

1 n n r 0 

n=o 

(7.8) 


N-l 


where 


T^ = total time interval of interest. 


S n 


y + 1 


2nn 

% 


W = e i2 * /N . 
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L 

B(n) = J Ln(f(x.r + i(n + IN) I s )} 

1-0 (7.9) 

Thus, the numerical values of f(x,T) are computed at N equally 
spaced time points = jAT , j = 0,1 ,2 ,...,N-l . For best results, it 
is suggested that the product LxN must range from 50 to 5000 and 
from 5 to 10 . The computations involved in equation (7.8) are performed 
by employing the Fast Fourier algorithm of Cooley and Tukey (1965). In 
case of the Fourier transform (s = — iu») , the above algorithm is equivalent 
to a Fourier synthesis. 

The above algorithm was tested for a number of trial functions (Ref. 
Ahmad, 1983). For L = 1, N ■ 200, and = 6 the numerical inversion 

results were highly accurate. Using N = 20 and neglecting the results for 
very early time steps (up to t = 0.05T) and for late time (after t = 0.75T) 
introduces a maximum error of only 2-3 percent and an average error of 0.6 
percent. Since use of N = 20 results in very substantial savings in 
computation time, this option is employed for three-dimensional problems 
and the results are plotted up to 15 time steps (i.e. T = 0.7 5T^). 
However, for two-dimensional problems both N = 20 and N = 50 are used. 


VII. 5 


In order to demonstrate the range and accuracy of the transformed 


domain solution with the numerical inverse transformation, a series of 


examples are presented ranging from a simply supported beam to a cavity in 
infinite space. The accuracy of the technique developed is compared to the 
available analytical and numerical results. In all cases, English units 
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are used with foot (ft) for length, pound (lbf.) for force, and second (s) 
for time, except otherwise specified. 

(A) Two-dimensional Applications 

(a) S iittplv-supoo rted beam subjected to step loading : 

A simply supported beam with a rectangular cross-section is subjected 
to a uniformly distributed step pressure as shown in figure 7.1. The 
dimensions of the beam are, length L = 30, depth d = 2, and width w = 1. 

7 

The material properties are, modulus of elasticity E * 3 x 10 , Poisson's 

—3 

ratio v — 0.3, and mass density p = 0.733 x 10 . The purpose of this 

analysis is to compare the solution predicted by the present method with 
that reported by Bathe et al (1974) by using NCNSAP. The Boundary element 
mesh as well as the finite element mesh are also shown in figure 7.1. 

Figure 7.2 shows the response (i.e. deflection at midspan) calculated 
using BEM and that from NONSAP. The time step used in the finite element 
solution to obtain the same results from the Wilson O and the Newmark 
integration schemes was AT = 0.5 x 10 -4 sec; whereas, the time step used 

_3 

in the present analysis is AT = 0.5 x 10 sec. In spite of the larger 
time step, the present analysis produces results identical to that reported 
in NONSAP. This help confirm the high accuracy and stability of the method 
presented in this chapter. 

(b) Half-space under prescribed time-dependent stress distribution : 

In this application, the results obtained by the present transformed 

domain, transient, dynamic formulation are compared against the solutions 
from finite difference by Tseng et al (197 5) and those from time-domain 
Boundary elements by Mansur and Brebbia (1985). 

The problem to be analyzed is depicted in figure 7.3(a). The half- 
space was initially at rest and then a part of its surface is disturbed by 
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a vertical pressure which is continuous in both time and space. Tseng used 
a transmitting boundary along with a generalized lumped parameter model to 
analyze this problem. His finite difference grid is shown in figure 
7.4(a). The boundary element discretization is shown in figure 7.4(b). 

The material properties of the half-space are, modulus of elasticity 
E = 200 ksi, Poisson's ratio v = 0.15 , and mass density 

« A J 

p = 1.9534 x 10“* lb-sec / in . For this problem, the time increments used 
by Tseng and Mansur and Brebbia was AT = 1 msec and AT = 3.65 msec , 
respectively whereas, in the present analysis, a much larger time 
increment, AT - 6 msec . is used. 

The time history of the vertical displacements plotted in figures 
7.3(b), 7.5, 7.6 and 7.7 are in reasonably good agreement with the previous 
results, even though a larger time- increment is used in the present 
analysis. The major difference in the results are in the displacements of 
point G(150,-10). in Tseng's work, this point is located on the 
transmitting boundary hence the finite-difference displacements at this 
point are not accurate. Similarly, in the case of boundary element 
analysis by Mansur, this point is located just belcw a boundary node which 
is a very difficult point to calculate interior displacements In the 
present analysis, none of the above mentioned problem is present and thus 
the displacements obtained in the present work is more accurate. The 
difference between the displacements, at point F( 80,-60) obtained by 
Mansur and Brebbia and present analysis is probably caused by the error due 
to numerical integrations. The present analysis uses a more sophisticated 
integration scheme than that used by Mansur and Brebbia and hence the 
results obtained by the present analysis should be more accurate. 

The time history of stresses at points A(45,-75), B(75,-75) and 
C(5,75) are plotted in figures 7.8, 7.9 and 7.10, respectively. It can be 
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seen that the results from the present analysis are in very good agreement 
with the results reported by Tseng during earler times, and are in good 
agreement with those reported by Mansur and Brebbia during later times. 
The difference at short times is due to an approximation used by Mansur and 
Brebbia in the calculation of interior stresses, i.e., the stress at a 
interior point is obtained by calculating the stresses on a triangular cell 
with the specified point as its centroid, whereas, the difference at later 
times is caused by errors generated at the transmitting boundaries used by 
Tseng. Finally, it should be noted that the results from the present 
analysis are in reasonably good agreement with the finite difference and 
the time-domain, boundary element solutions. 

(c) Semi-infinite beam subjected to a suddenly-applied 
bending moment : 

A semi-infinite beam simply supported along its edge is subjected to a 
suddenly applied bending moment M = m q H(T- 0), as shown in figure 7.11. 
The beam is considered to be under a plane-stress condition and the Poisson 
ratio is taken as v = 1/3 . 

A finite element analysis of this problem was carried out by Fu 
(1970), and a boundary element analysis was carried out by Mansur and 
Brebbia (1985). Boley and Chao (1958) obtained the results for the same 
problem using beam theory. Transverse displacements along the axes of the 
beam obtained by the above researchers and the present method are shown in 
figure 7.12. This displacements plotted in figure 7.12 refer to T = 5r/c Q 
where r is the radius of gyration of the beam cross section and c Q is 
the one-dimensional wave propagation speed. 

In the present analysis two types of boundary conditions are used. In 
the first case, the beam is fixed from transverse movement by incorporating 
zero transverse displacement at the midpoint of the finite end of the beam. 
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and in the second case, zero transverse displacements are incorporated for 
all the nodes along the finite end of the beam. The displacements obtained 
by incorporating the first boundary-condition case are in good agreement 
with finite element results whereas the displacements obtained by using the 
second boundary-condition case are in good agreement with the beam theory 
and Mansur's solutions . Therefore, the difference in the results of 
Mansur and Fu are essentially due to end boundary conditions. 

(B) Three-dimensional Applications 

(a) Cantilever Beam. subjected to time-harmonic axial, tension : 

A uniform beam with a rectangular cross-section has a modulus of 
elasticity E = 1.16 x 10 7 , a Poisson's ratio v = o.o , and a mass 
density p = 2.0 . It is fixed at one end and a uniformly distributed 
axial tension p = 1000 sinQT , Q » 0.628 r/s , is applied at the free end. 
Traction-free conditions hold along the sides. The dimensions of the beam 
are length L = 4 , depth d = 2 , and width w = 1 . The beam is modelled 
by six quadrilateral elements resulting in 20 nodes, as shown in Fig. 7.13. 
The same figure plots the axial displacement at the free end as a function 
of time along with the analytic solution developed from Clough and Penzien 
(1975). Agreement is very good considering that only 20 points were used 
in the Laplace transform domain and that the sinusoidal load was 
represented by straight line segments for the purpose of the direct Laplace 
transformation, Eq. (7.5). 

(b) Spherical cavity in infinite space: 

A spherical cavity is embedded in an infinitely extending medium with 
E = 8.993x10^ , v = 0.25 , and p = 2.5xl0 -4 . The radius of the cavity is 
a = 212 and its surface is discretized into 3 triangular elements per 
octant for a total of 50 nodes, as shown in figure (7.14). The 
characteristic times required for the pressure and shear waves to travel a 
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cavity radius are 0.00102s and 0.00177s , respectively. Two cases are 
considered : 

(i) Spherical cavity under sudden radial pressure : A radial pressure 
p = 1000 is suddenly applied and maintained at the cavity surface. Figure 
7.15 shows the radial displacement history obtained by using the inverse 
Laplace transform algorithm with only 20 data points. The response is 
obtained for two different time steps, AT equal to 0.0005s and 0.0003 5s . 
Concurrently plotted is the exact solution (Ref. Timoshenko, 1970). In 
general, the numerical results are in good agreement with the analytical 
solution. These is some oscillation in the Laplace transform solution 
towards the end of the total time so that about 85% of the time spectrum 
obtained is actually plotted. 

(ii) Spherical cavity, engulfed bv. a pressure wave : A propagating 

plane pressure wave whose front is perpendicular to the Z-axis first 

impinges on the pole with coordinates (0,0,212). The resulting non-zero 

incident stresses are “ -1000 H(T-T q ) , and = 

(W(l-\>))<r , where H is the Heaviside function and T the time 
zz o 

required for the wave to reach the station in question. This wave 
propagation type of problem is solved by superposition (Ref. Eringen, 
1975). A three-quadrilateral s-per-octant mesh resulting in 74 nodes (Ref. 
figure 7.14) is used here in conjunction with the numerical inverse 

transformation utilizing 20 data points. Figure 7.16 shows the hoop 

* * 

stresses , and normalized by the magnitude of the incident stress 

a zz* } versus the non-dimensional time x* = aT/c 1 . The plots are for 
three locations on the surface of the cavity: the two poles (6 = 0,jr) and 
the equator (d = n/2) . Concurrently plotted are the analytic results 
(Ref. Pao and Mow, 1973). Good agreement is observed between the two 
solutions. 
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VII.6 CONCLUDING REMARKS 


An advanced algorithm based on the transformed domain boundary element 
formulation for transient dynamic analysis has been presented. The 
numerical implementation employed is one of the most general currently 
available and can be used in conjunction with substructuring to treat two 
and three-dimensional solids of complicated geometry and connectivity. 
Interior, exterior and halfspace problems can all be solved by the present 
algorithm. The current implementation is also capable of handling sliding 
interfaces in the soil-structure interaction problems. Thus, the algorithm 
presented is a viable alternative to that based on finite element 
methodology . 
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CHAPT ER. V-III 

TIME DOMAIN TRANSIENT DYNAMIC ANALYSIS 
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VIII. 1 INTRODUCTION 


The work described in this chapter is based on the numerical 
implementation of the direct boundary element method for time-domain, 
transient analysis of three-dimensional solids in a most general and 
complete manner. The present formulation employs the space and time 
dependent fundamental solution (Stoke's solution) and the Graffi's cfynamic 
reciprocal theorem to formulate the boundary integral equations in the time 
domain. A time-stepping scheme is then used to solve the boundary- initial 
value problem by marching forward in time. Interpolation functions in 
space and time are used to approximate the field quantities, and a 
combination of analytical (time-integration) and numerical integration is 
then carried out to form a system of linear equations. At the end of each 
time step, these equations are solved to obtain the unknown field 
quantities at that time. 

In the following sections, a description of the proposed methodology 
is presented in detail. The materials related to the representation of 
geometry, spatial variation of field quantities, numerical integration and 
solution of equations at each time step are similar to those already 
described in Chapter VI for one value of transform parameter s except for 
the fact that, in the present case, all the quantities are real. The 
matrix equation solver used for the present case is a real-variable version 
of the out of core complex solver described in Chapter IV, Sec. 4.G. The 
built-in symmetry and substructuring capabilities described in Sec. VI.3.B 
are also included in this implementation. A number of numerical examples 
are finally presented to demonstrate the stability and accuracy of this 
dynamic analysis technique. 
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VIII. 2 TRANSIENT BOUNDARY INTEGRAL FORMULATION 

The direct boundary integral formulation for a general, transient, 
elastodynamic problem can be constructed by combining the fundamental 
point-force solution of the governing equations (4.1) (Stoke's solution) 
with Graff i's dynamic reciprocal theorem. Details of this construction can 
be found in Banerjee and Butterfield (1981). For zero initial conditions 
and zero body forces, the boundary integral formulation for transient 
elastodynamics reduces to: 

C ij ( £ )u i (i ’ T) = I [G ij ( X.£.T)*t.(x,T) 

S 

“ F i j(x.£.T)*u i (x.T)ldS(x) (8.1) 

where: 

r T 

G ij* fc i = J G ij(x»T;£,T)t^(x.T)dT 

o 

r T 

F^j*u^ = J F^j (x,T;£,T)u^(x,T)dT (8.2) 

O 

are Riemann convolution integrals and £ and x ate the space positions of 
the receiver (field point) and the source (source point). The fundamental 
solutions and F^j are the displacements and tractions at a point x 

and at a time T due to a unit force vector acting at a point L at a 
time r . These functions are listed in a compact form in Appendix A4. 

Equation (8.1) represents an exact formulation involving integration 
over the surface as well as the time history. It should also be noted that 
this is an implicit time-domain formulation because the response at time T 
is calculated by taking into account the history of surface tractions and 
displacements up to and including the time T . Furthermore, equation 
(8.1) is valid for both regular and unbounded domains. 
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Once the boundary solution is obtained, the stresses at the boundary 
nodes can be calculated without any integration by using the scheme 
described in Sec. VI.3.D. For calculating displacements at interior points 
equation (8.1) can be used with c^ = 5^ and the interior stresses can 
be obtained from 

<T jk (£,T) - J tG? jk (x.I,T)*t i (x.T) - Fjj k (x.i.T)*U i (*,T)]dS(x> 

S (8.3) 

The functions G® j k and F^j k of the above equation are listed in 
Appendix A5. 

The constitutive equation and the boundary and initial conditions are 
described in Chapter IV (Ref. Secs. IV. 1 and IV. 2). 


VIII. 3 TIME STEPPING SCHEME 

In order to obtain the transient response at a time T^ , the time 
axis is discretized into N equal time intervals, i.e. 

N 

T^ = J nAT (8.4) 

n=l 


where at is the time step. 

Utilizing equations (8.4) and (8.2), equation (8.1) can be written as: 


,N 


c ijV^V - 1 f ‘Vi 


Vi s 


F^^ldSdr 




[G. -t. - F. .u. IdSdt 

l] l l] l 


(8.5) 


t=o 
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where the integral on the right hand side is the contribution due to past 
dynamic history. 

It is of interest that equation (8.5) like equation (8.1) still 
remains an exact formulation of the problem since no approximation has yet 
been introduced. However, in order to solve equation (8.5), one has to 
approximate the time variation of the field quantities in addition to the 
usual approximation of spatial variation. For this purpose two types of 
interpolation functions are used which are described with the resulting 
time-stepping algorithms as follows. 

(A) Constant. Time Interpolation 

In this case, both displacements and tractions are assumed to remain 
constant during a time step, i.e., 

N 

u ^(x.t) * £ u!?(ji)<6 n (T) and 
n=l 

N 

VX.*) = J tj(x>d n (r) (8.6) 

n=l 

where 

^ n (r) = 1 for (n-l)AT < t < nAT , and 
= 0 otherwise; and 

u?(x) and t?(&) represents the spatial variation of u^ and t^ , 
respectively, at time T n . 

For illustrative purposes, first consider the form of equation (8.5) 
for the first time step; i.e. 

T i 

C ij u i'i- T l> - J J Wyti - r^Mafe - 0 (8.7) 

T s 
o 
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The time integration in (8.7) is done analytically (Ref. Appendix Dl) and 
the surface integration is performed in the usual manner (i.e. 
numerically). After the integrations and the usual assembly process, the 
resulting system of algebraic equations is of the form: 

[A 1 ] tX 1 } = [B 1 ] {Y 1 } (8.8) 

where A and B are coefficient matrices, Y and X are the known and 
unknown components of the boundary tractions and displacements 
respectively, and the superscript pertains to the time step. 

New consider the boundary integral equation for the second time step, 

i.e. ; 

T 2 

C ij“i'*- T 2> - J J - P ij u i ! 

Tj S 

T 1 

= J | iG^ti - F ijUi ]dSdT (8.9) 

T q S 

If the time interval (T 2 ~t 1 ) is same as (Tj-T 0 ) the resulting 
coefficient matrices of the left hand sides of equations (8.7) and (8.9) 
become identical. This is so because the time translation properties of 
the fundamental solutions and (Ref. figure 8.1), contain time 

functions with arguments (T— c) and therefore the convoluted integral 
corresponding to the interval T x <_ x £ T 2 with T = T 2 is identical to 
that of the interval T q <_ x <. T with T = T t . 

The right hand side of equation (8.9) is evaluated at time T = T 2 
with the time integration over the interval T Q to T 1 and thus provides 
the effects of the dynamic history of the first time interval on the 
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current time node (i.e. T 2 ). mow, the resulting system equation for this 
time node (T 2 ) is of the form: 

[A 1 ] tX 2 } = [B 1 ] {Y 2 } - [A 2 ] {X 1 } + [B 2 ] CY 1 } (8.10) 

where superscripts on X and Y pertains to the time nodes and superscripts 
on A and B denote the time step in which they are calculated. 

Using equations (8.8) and (8.10), equation (8.5) can be written in an 
assembled form as: 

N 

[A 1 ] {X^} = [B 1 ] (Y N ) - ]> [[A n HX N_n+1 ] - [B n HY N ' n+1 )] 

n=2 (8.11a) 

or [A 1 ]^) = [B 1 ] (Y N ) + (R N J (8.11b) 

where R N is the effect of the past dynamic history on the current time 
node . 

The above equation can be solved to find the unknown X N at time T^. 
It may appear at first glance that a prodigious coefficient calculations 
are involved. However, a closer examination will reveal that: 

(i) If the time step size is constant, the A 1 and B 1 matrices do 
not change from time step to time step. 

(ii) For each time step, a new R N needs to be formed. This 
involves the evaluation of a new set of coefficients A n and B n involving 
the effects of the dynamic history of the first time interval on the 
current time node. Eventually, however, this contribution to R N reduces 
to zero and from that point onwards no new coefficients need to be 
evaluated . 
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In the present implementation, the representative values of the 
displacements and tractions during a time stepping interval is obtained by 
averaging the values of these quantities at two time nodes of that 
interval . 

(B) Linear Time Interpolation 

In this case, both displacements and tractions are assumed to vary 
linearly during a time step.i.e. 


N 

U^x.t) = J f^-V) +H 2 uJ( X )] 
n=l 

N 

J [M 1 tJ _1 (x) + M 2 tJ(x)] (8.12) 

n=l 


where and M 2 are the time functions, and are of the form: 


T - r 

M = ~ 6 (t) 

l AT n 


M 


2 



V r) 


(8.13) 


Again for illustration proposes, consider the boundary integral equation 
for the first time step, i.e. 

T i 

c ij u i (5 ' T l ) ‘ ] I [G ij fc i _ Fj.j^ldSdr = 0 (8.14) 

T o S 

The time integration in equation (8.14) by utilizing (8.12) is done 
analytically (Ref. Appendix D2). After the usual numerical integration and 
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assembly process, the resulting system equation is of the form: 

tA^nx 1 } -[B^HY 1 } + [aJhX 0 } - [bJhy°} = {0} (8.15a) 

where: 

A and B are the matrices related to the unknown and known field 
quantities, respectively; 

X and Y are the vectors of unknown and known field quantities, 
respectively; 

for X and Y superscript denotes the time; 

for A and B supercript denotes the time step at which they are 
calculated, and the subscript denotes the local time nodes (1 or 
2) during that time-stepping interval. 

Since all the unknowns at time T = 0 are assumed to be zero, 
equation (8.15a) reduces to 

[A^HX 1 ) - IB^HY 1 } + [bJhy 0 } (8.15b) 

For second time step, the assembled system equation has the form 
[A*](X 2 ) - [B 2 HY 2 } + [aJhx 1 } - fB^HY 1 } = 

- HA^HX 1 ) - [B^HY 1 } + [A 2 HX°) - [B 2 HY°}} (8.16a) 

Similar to the constant time variation scheme, only the matrices on the 
right hand side of equation (8.16a) need to be evaluated. However, one 
needs to integrate and assemble four matrices at each time step as compared 
to two in the case of constant time variation. This can be done with only 
a small increase in computational time by integrating all the kernels 
together and then assembling all the matrices together. Equation (8.16a) 
can be rearranged such that: 
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(8.16b) 


[A^nx 2 ] = [B 2 HY 2 } - [aJ+A^HX 1 } - [bJ+B^HY 1 } + [B 2 HY°} 

In the above equation, all the quantities on the right hand side are known. 
Therefore the unknown vector X at time T 2 can be obtained by solving 
the above equation. 

Thus, for the present case, the boundary integral equation (8.5) can 
be written in a discretized form as: 

N 

(A 2 nx N ) - [b 2 hy n ) - - J [[aJ+aJ - 1 ] (x N_n+1 } 

n=2 

- [B5+Bj" 1 HY N “ n+1 }] + [bJHY 0 ) (8.17) 

or 

IAJ] fX N ) = [B^lfY 11 } + (R N ) (8.18) 

The discussion in the previous section regarding the causal properties 
of the fundamental solution holds true for the present case also. 

It is of interest to note that, if time interpolation functions Mj 
and M 2 are replaced by = M 2 = 0.5 d n (r) , the time stepping scheme 
for linear variation can be used for the case of constant variation with 
averaging between the local time nodes. 

VIII. 4 SOME ASPECTS OF NUMERICAL IMPLEMENTATION 

The numerical implementation of the boundary integral equation for 
time-domain, transient elastodynamics is essentially similar to that 
described in Chapter VI for steady-state elastodynamics, except for the 
following : 

(i) All the quantities involved in the time domain analysis are real, 
instead of complex as in the case of steady-state dynamics. 
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(ii) There is a fundamental difference between static or steady-state 
dynamic analysis and time-domain transient analysis when it comes to the 
numerical integration schemes outlined in Chapter VI. In the static or the 
steady-state case, the integrands in all of the nonsingular surface 
integrals are infinitely differentiable and solution accuracy can, 
therefore, always be improved by the use of increased integration order. 
In the transient case, however, the point load solutions are only 
continuous. Physically this corresponds to the fact that the disturbance 
at some later time due to an impulse applied to a spatial location at a 
given time (past) is only present in a finite portion of the space (Ref. 
figure 8.19). This means that the kernel function may be nonzero over only 
part of a given surface element. While the integrand is infinitely 
differentiable within both the zero and nonzero regions considered 
separately, its overall continuity over the entire element is only C . 

i w 

The use of higher order quadrature rules is, therefore, of little use in 
improving accuracy. Based on these observations, a revised integration 
strategy was adopted for the transient case. All surface elements are 
subdivided into a relatively large number of subelements and relatively 
low-order (usually 2nd or 3rd) quadrature rules together with the usual 
distortion in mapping (so that the kernel shape functions and Jacobian 
products remain well behaved) are used over each subelement. This has led 
to much improved accuracy in the transient analysis. 

(iii) In the case of singular integration, the subelements are 
subdivided in the radial direction also. This subdivision has been found 
to increase both the accuracy and the stability in the time domain 
approach. 

(iv) In time domain analysis, the fundamental solution as well as the 
field variables are functions of real time T and therefore the system 
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equation at each step corresponds to a time T rather than to a 
transformed parameter s as in the case of Laplace domain analysis (Ref. 
Chapter VII). 

(v) All the matrices related to the past convolution are stored on 
sequential tapes, and at each time incr orient they are used along with the 
boundary excitation history (of tractions and displacements) to calculate 
the effect of the past dynamic history on the current time node. 

(vi) The matrix B 1 (Ref. eq. 8.11) is stored on a sequential tape arid 
at each time increment it is used to calculate the contribution to the 
right hand side due to the known field quantities at the current time. 

(vii) During the solution process at the first time increment, the 
decomposed form of the system matrix A 1 (Ref. eq. 8.11) is stored on a 
direct-access file for later use. After that, at each time step, all of 
the known tractions and displacements are multiplied by appropriate 
coefficient matrices to form a new right-hand-side vector. The decomposed 
form of the system natrix is then used with the new right hand side vector 
to calculate the unknown displacements and tractions at the current time. 
This process of repeated solution by using the decomposed form of the 
system matrix is highly efficient and thus results in considerable saving 
in solution time. 


VIII. 5 


NUMERICAL ACCURACY. STABILITY AMD CONV 




:e of solution 


In order to investigate the accuracy, stability and convergence of the 

proposed numerical technique, the problem of the radial expansion of a 

spherical cavity in an infinitely extending medium, subjected to suddenly 

applied and maintained internal pressure [p(T) = 1000] was studied. The 

material properties were as follows: E = 8.993xi0 6 psf, v = 0.25 , and p 

—9 2 4 

= 2.5x10 lb-sec /ft . The radius of the cavity was taken as R = 212 ft 
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and three different meshes shown in figure 7.11 were used to discretize the 
cavity surface. Using built-in symmetry capabilities, this problem was 
modeled by one octant only. The first mesh has one six-noded triangular 
element, the second has three triangular elements (total of 10 nodes), and 
the third has three eight-noded quadrilateral elements (total of 16 nodes). 

In figures 8.2-S.4 the radial displacement u f (r = R.T) normalized 
by the static value is plotted against time for a total of nine different 
time steps. Concurrently plotted is the exact solution (Ref. Timoshenko, 
1970). These results conclusively demonstrated the unconditional 
stability of the BEM formulation. The accuracy is highest when the time 
step is between 1/3 to 3/4 of the characteristic time R/Cj . In all 
cases the results approach the static response without exhibiting any 
supurious oscillations. 

The effect of the surface discretization is demonstrated in figure 
8.5, where the time variation of u r (R»T) is plotted for all three meshes 
for the same time step AT = 0.0003 5 s . It is observed that the errors 
in the dynamic response are consistent with the average error committed in 
the static response which is 12% for the first mesh, 3% for the second 
mesh and 1.5% for the third mesh. Thus, the numerical technique presented 
here converges to the actual results with finer discretization of the 
surface of the boundary. 

VIII. 6 EXAMPLES OF APPLICATIONS 

A number of representative problems are chosen to test the accuracy 
and the stability of the time-stepping solution. In all cases, English 
units are used with foot (ft) for length, pound (lbf) for force, and 
second (s) for time. 
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(a) Bar Subjected to a Transient End Load . 

(i) Square cross-section : A bar with square cross-section is held 

along its sides by lubricated rollers and is fixed at one end. The free 

end is subjected to a suddenly applied and maintained uniform compression 

t = 1000. The dimensions of the bar are L = 8.0 and b = 2.0 . in view 
z 

of the material properties, the characteristic time required for the 

compressive wave to reach the fixed end is 0.03 578 sec. Figure 8.6 shows 

the discretization and the numerical results for the normal stress a . 

zz 

in which the results from the time domain algorithm for two different time 
steps AT are compared with the exact analytical solution for one- 
dimensional stress wave propagation (Ref. Timoshenko, 1970). Although the 
numerical results are in good agreement with the analytical solution, it is 
clearly very difficult to reproduce the sharp jump in the stress as the 
disturbance reaches the point initially and when the reflected stress wave 
returns to the same location. This difficulty has been observed elsewhere 
as well (Ref. Belytschko et al, 1976). 

The axial displacement history at the free end is shown in figure 8.7. 
The displacements are normalized by static displacements and the time is 
normalized w.r.t the characteristic time required for the compressive wave 
to reach the fixed end. It can be seen that the numerical results are in 
good agreement with the analytical solution. The differences are mainly 
due to the three-dimensional nature of the simulated problem. 

(ii) Circular-cross-section : In order to investigate the effects of 

the cross-section on the numerical results, a bar with circular cross- 
section having the same material properties and boundary conditions as 
described in the last example was analyzed. The boundary element mesh for 
this problem is shown in figure 8.8. The bar has a length L = 5 and 
diamater d = 1 . Thus, the characteristic time required for the 
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compressive wave to reach the fixed end is 0.02236 sec. The time step used 

in this example is AT = 0.004475 sec. 

Figure 8.9 shows the numerical results for the normal stress a at 

zz 

the midspan of the bar against a one-dimensional analytical solution. As 
mentioned in the last example, the sharp jumps in stress are diffused in 
the numerical results. However, by using more elements and smaller time 
steps, the numerical results in the vicinity of the jumps will agree more 
closely with the analytical solution. 

The time history of the normalized axial displacements at the free-end 
is plotted in figure 8.10 against the one-dimensional analytical solution. 
The results are in good agreement, except for the peak displacements. The 
numerical peak values are less than that of the analytical solution and 
this results in an increase in the difference between the two solutions at 
later times. The difference, once again, is mainly due to the three- 
dimensional nature of the problem under consideration. 

(b) Spherical Cavltv. . 

A spherical cavity is embedded in an infinitely extending medium with 
E = 8.993xl0 6 , v = o.25 , and p = 2.5xl0 -4 . The radius of the cavity R = 
212 and three different meshes for its surface discretization are shown in 
figure 7.11. Using the built in symmetry capabilities, this problem is 
modeled by one octant only. The characteristic times required for the 
pressure and shear waves to travel a cavity radius are 0. 00102s and 
0.00177s , respectively. Four cases are considered: 

(i) Spherical cavity under sudden radial expansion : A radial 

pressure p = 1000 is suddenly applied and maintained at the cavity 
surface. Figure 8.11 shows the time variation of deviatoric stress at the 
cavity surface obtained by the time domain algorithm. Concurrently plotted 
is the result reported by Hopkins (1960) based on the work of Hunter 
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(1954). In general, the numerical results are in good agreement with the 
analytical solution. The transient, time-domain solution remains stable 
and reaches the expected static solution at larger times. It can be ssen 
that the maximum deviatoric stress for the transient case is 1.77 times the 
applied pressure whereas for static case it is 1.54 times the applied 
pressure. 

(ii) Spherical cavity subjected to a rectanaualr pulse of radial 
pressure : A triangular pulse of radial pressure, as shown in figure 8.12 
is applied at the cavity surface. This example is solved by using linear 
time interpolation functions and two different time steps. The radial 
displacements at the cavity surface are plotted in figure 8.12. The 
numerical results from both the time steps are almost identical. Thus, 
this example once again demonstrates the stability of the present 
algorithm. 

(iii) Spherical cavity subjected to a rectangular pulse of radial 
pressure : A rectangular pulse of radial pressure as shown in figure 8.13 
is applied at the cavity surface. This example is also solved by using 
linear-time interpolation functions and two different time-increments. 
Figure 8.13 shows the time history of the radial displacement of the 
cavity. By comparing this results with those due to a triangular pulse 
(i.e. fig. 8.12), it can be seen that, in general, displacements at any 
time interval due to the rectangular pulse are twice that due to the 
triangular pulse. This is because the response depends upon the total 
inpulse and the total inpulse due to the rectangular pulse is double that 
due to the triangular one. Hence, the displacement amplitude response due 
to the rectangular pulse is also approximately double that of the 
triangular pulse. In addition, since the energy input is the same in both 
problems, the response curves for both cases have the same shape. 
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(iv) Spherical ca v ity engulfed bv a pressure wave : A propagating 

plane pressure wave whose front is perpendicular to the Z-axis first 

impinges on the pole with coordinates (0,0,212). The resulting non-zero 

incident stresses are a ^ = -1000 , <s ^ = a ^ = (v/(l-v))<r 

zz xx yy zz 

This wave propagation type of problem has been solved by superposition 
(Miklcwitz, 1978). The three quadrilaterals per octant mesh is used here 

in conjunction with the time-domain approach. Figure 8.14 plots the hoop 

* * 

stresses and <x eo normalized by the magnitude of the incident stress 

a ^ versus the non-dimensional time -c* = RT/c. . The plots are for 
three locations on the surface of the cavity: the two poles (6 = o,n) and 
the equator ( 6 - nil) . Concurrently plotted are the analytic results (Pao 
and Mow, 1973; Norwood and Miklowitz, 1967), obtained by analytical 
inversion of the Fourier transformed solution. Good agreement is observed 
between the two solutions. Finally, figure 8.15 plots the radial 
displacement time history at the same three locations as before. 

(c) Transient point load on halfsoace . 

This example is Lamb's problem for an impulsive vertical point force 
on the surface of a semi-infinite solid (Pekeris, 1955). The modelling 
difficulty encountered here is that the point load must be represented by a 
finite (yet small) area, and hence the complicated mesh shown in figure 
8.16. There are 7 co-centric rings, with four elements per ring, resulting 
in a total of 85 nodes. However, by using the symmetry only one quarter of 
the mesh is modeled as a input geometry data. Uniform vertical tractions 
t = looo are prescribed on the triangular elements of the inner ring 

Z 

which has an outer radius of 0.0 5. The outer ring has an inner radius of 
3.0 and is modelled by infinite elements. Obviously the small circular 
load solution behaves differently from the analytical, point-force solution 
(labelled solution A). The static solution showed that the results agree 
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well at a radial distance of 0.2 where there is a 4.5% and a 1.0% 
difference in the horizontal and vertical displacements, respectively. 
Therefore, the results from the time domain algorithm shown in figure 8.17 
are for the normalized horizontal displacements at r = 0.2 . when the 
exact solution (Pekeris. 1955) is used to calculate the superimposed 
effects of multiple point forces to reproduce a finite area loading 
(labelled as analytic solution B), good agreement with the BEM results is 
obtained. 

(d) Square flex ible foot ing on half-space . 

In this example, a square flexible footing on half-space is subjected 
to a time dependent vertical tractions. The mesh for this problem is shown 
in figure 6.8(a) and is the same as that used for calculating vertical 
compliance for rigid square footing. The side of the footing is B = 2b = 
2 , and the material properties of the half-space are: elastic modulus E 
= 2.6, Poisson's ratio » = 0.3 and mass density p = l.c . The time step 
used for this analysis is AT = 0.2 . The time history of the applied 
pressure and the vertical displacements at the center and corner of the 
loaded area are plotted in figure 8.18. It can be seen that the vertical 
displacement at the center of the footing converges to the static value 
after 2.4 seconds. Whereas the vertical displacement at the corner of the 
footing seems to be converging to its static value at a later time. The 
mesh used for this problem gives a maximum error of 2% for static 
analysis, hence the results obtained for the present problem are supposed 
to be reasonably accurate. Finally, this example shows the usefulness of 
the present algorithm for transient dynamic analysis of half-space 
problems. 
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VIII. 7 CONCLUDING REMARKS 

An advanced algorithm based on the direct boundary elanent formulation 
for time-dependent elastodynamic analysis of three-dimensional solids has 
been presented. The algorithm is an unconditional ly-stable, implicit, 
time-marching scheme and is capable of producing very accurate results. 
However, for better accuracy, it is recommended that the time step should 
remain smaller than L/Cj , L being the smallest distance measured along 
the surface between two corner nodes of an element. This algorithm is a 
viable alternative to that based on the finite element methodology, 
particularly for soil-structure interaction problems. 
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CHAPTER IX 

NONLINEAR TRANSIENT DYNAMIC ANALYSIS 
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IX. 1 INTRODUCTION 


In this chapter, a direct boundary element formulation and its 
numerical implementation for nonlinear transient dynamic analysis of 
three-dimensional deformable solids of arbitrary shape and connectivity is 
presented. The formulation is based on an initial stress approach, and is 
the first its type in the field of Boundary Element technique. The 
nonlinearity considered in this analysis is that due to the nonlinear 
constitutive relations, i.e. material nonlinearity. The boundary integral 
equations are cast in an incremental form, and thus, elasto-plastic 
relations of the incremental type are used for material description. These 
equations are solved by using a time-stepping algorithm in conjunction with 
a iterative solution scheme to satisfy the constitutive relations. The 
resulting algorithm is an unconditionally stable implicit scheme. However, 
the size of the time step that can be used is restricted by the size of the 
elements used for modelling the surface of the problem under consideration. 

In the present analysis, the geometry and the field variables are 
represented by higher-order isoparametric shape functions to model complex 
geometries and rapid functional variations accurately. In this chapter, 
the discussion first focuses on the formulation of the method, followed by 
the numerical technique for discretization and spatial integration of 
volume integrals. For discretization and spatial integration of surface 
integrals, the numerical integration techniques developed in earlier 
chapters (Ref. Secs. VI.3 and VIII.4) are used. The material pertaining to 
the time-stepping scheme along with the iterative solution algorithm are 
presented next. Numerical examples are finally presented to demonstrate 
the accuracy and applicability of the present method. 
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IX. 2 BOUNDARY INTEGRAL FORMULATION FOR DYNAMIC PLASTICITY 

The direct boundary integral formulation for a nonlinear transient 
dynamic problem, based on an initial stress approach, can be constructed by 
following a procedure similar to the one that has been used for a nonlinear 
static problem (Ref. Sec. 12.4(b), Banerjee and Butterfield, 1981). Under 
zero initial conditions and zero body forces, the boundary integral 
equation for nonlinear transient dynamics is of the form 

CijdJUid.T) = J [Gy (x,£.,T)*t i (x,T) - F i j(x,I.T)*u 1 (x.T)]dS(x> 

S 

+ j B il j(x,5.,T)*o? 1 (x.T)dV(2i) (9.1) 

V 

where * denotes convolution (Ref. Sec. VIII. 2); 

i and x are the space positions of the receiver (field point) and 
the source (source point), respectively; 
is the initial stress tensor; 

V denotes the volume of the body; and 

the fundamental solutions Gy, and By^ are listed in 

Appendices A4 and A 6 . 

Assuming all the field quantities to have a zero value at time T = 0, 
the boundary integral equation (9.1) can be written in an incremental form 
as follows: 

Cy^Au^I.T) - | [Gy(x,l,T)*At i (x.T) - Fy (x,£.,T)*Au i (x,T) IdS(x) 

S 

+ I B iij(X»I.T)*A0? 1 (ji.T)dV(x) (9.2) 

V 

where A denotes the incremental quantity. 

The stress increment at an interior point £. can be obtained by 
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taking derivatives of equation (9.2) and using the constitutive 
relationships (Ao„ = }<l Ae kl " Aa ij 5 aS: 

A«j k (I.T) = J [G? jk (x,I,T)*At i (x.T) - F®j, < (x,I,T)*Au i (x»T)]dS 
S 

+ J B iljk (x '£' T) * Ac il (x ' T) + J iljk Aaf 1 (l.T) 

V (9.3) 

The functions G?j k , F ijk' B iljk J iljk are Appendices A5 

and A 6. 

In equation (9.3), the volume integral must be evaluated in the sense 
of (V - V g ) with limit V g -> o and the tensor J.^ k is the jump term 
arriving from the analytical treatment of the integral over V g . This 
jump term is the same as that of static plasticity and is independent of 
the size of the exclusion V g provided the initial stress distribution is 
locally homogeneous (Ref. Banerjee and Davies, 1984; Raveendra, 1984; 
Banerjee and Raveendra, 1985). 

The equations for incremental stresses cannot be constructed at the 
boundary points by taking the field point (JL) in equation (9.3) to the 
surface due to the strongly singular nature of the integrals involved. 
However, the equations for incremental stresses at boundary points can be 
constructed by using a scheme similar to that described in Sec. VI.3.D. 
Using this scheme, the incremental stresses and the global derivatives of 
the incremental displacements at a boundary point JL b can be obtained by 
coupling the following set of equations: 

Aa i:j (i b .T) = [x.8 ij Au m>m (l b ,T) + n{Au Lj (I b ,T) + A Uj # . (I b ,T) } ] - Ac£.(£ b .T) 
At^.T) = A<Ty(£ b ,T)nj(£ b > , and 
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9Au . 


(9.4) 


l 


9ll a 


9AU. 

3^ 3l 1 a 


l_ 

where 11 is a set of local axes at the field point (1°) . 

The above equations can be combined together and written in a matrix 
form as 


[Slip) = tq) (9.5) 

where [SI is a 15x15 matrix which contains unit normals, a 3x3 unit 
matrix and material constants; p is the unknown vector of A<r.^ and 
9Au^/3?j ; and q is a vector containing the tractions At^ and local 
derivatives of the displacements Au^ . 

By making use of equation (6.5), the right hand side of equation (9.5) 
can be written as 

tq} = [Eltg] (9.6) 

where [E] is a 15x48 matrix of shape functions and derivatives of shape 
functions; and g is a vector of incremental nodal tractions and 
displacements over all of the local element nodes. 

Inverting matrix [S3 and utilizing equation (9.6), the set of 
equations (9.5) can be rearranged to form 

A« jk (i b .T) - 5? jk At l(i b .T) - ?J jk Aa l( l b .T) 

+ 5 iljk 4l ’il ( i b ' T) <9 - 7) 

It should be noted that the above equation is free of any integration and 
time convolution. 
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IX. 3 CONSTITUTIVE MODEL 


In dynamic plasticity, the choice of an appropriate constitutive model 
depends largely on the material properties and the loading conditions of 
the problem in hand. For this reason various constitutive models have been 
used for dynamic plasticity. However, for simplicity in the present 
analysis, the Von Mises model with isotropic variable hardening is used. 

In this model, the behavior in the elastic and plastic region is 
governed by the stress-strain relations: 


[ V 

Ae • • + . 1 . 5 . . 
n l-v n 


Ae. 


ifiAi 


« 1 - V « kk 


At u3 


(9.8) 


where A<r^ = D^j^Ae^ = incremental stress tensor, 


= incremental elastoplastic material modulus, 


4 = elastic shear modulus, 

<t q = uniaxial yield stress = V3S^S i j/2 , 

Sy = deviatoric stress = - 6 i j cr kk / 3 , and 

H = plastic-hardening modulus, the current slope of the uniaxial 
plastic stress-strain curve. 

The present implementation is suqh that any other constitutive model 
can be included without any difficulty. 


IX. 4 DISCRETIZATION AND SPATIAL INTEGRATION OF THE VOLUME INTEGRALS 
(A) Discretization 

Equations (9.2) and (9.3) provide the formal basis for developing the - 
dynamic plasticity algorithm. However, the initial stresses Act?^ defined 
in equations (9.2) and (9.3) are not known a priori and have to be 
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determined by satisfying the constitutive relations discussed in Section 
IX. 3. Thus, equations (9.2) and (9.3) and (9.8) can be regarded as a 
coupled system of nonlinear equations. In the present implementation, 
equation (9.3) and (9.7) are used to calculate the stresses at interior and 
boundary points and the nonlinear material model is then used to evaluate 
the inelastic stresses. Since the volume integrals of inelastic stress 
vanish except in regions of nonlinear material response, approximations of 
geometry and field quantities are required only where nonlinearity is 
expected. In the present work, isoparametric (quadratic) volume cells are 
used for approximating the geometry and the variation of initial stresses 
such that: 


x i M p* a)x ip 

«ij - (9 - 9) 

where x^ are carter ian coordinates, 

are nodal coordinates of the volume cell, 

Mp is a quadratic shape function for the volume cell, 
p represents the nodal points of the volume cell, and 
- denotes nodal quantities. 

A typical volume cell is shown in figure 9.1. 

The volume integral of equation (9.2) can be then represented as 

T 

I I B iii ? ( X.T;£ b ,r)A(j? 1 (x,r)dVdT 
o V 

L T 

= J J | B il j[s m (ll).T;I b ,r)]Mp( 1 i)dV m Acr?^ dx (9.10) 

m=l o V 
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where : 

y. 

Sr is the field point on the boundary (boundary node), 
a?*(2L) is the point in cell m , 

A«r?ip are the nodal values of incremental initial stress of the 
m*"* 1 cell, 

V m is the volume cell, and 

L is the total number of cell in a single region. 

Similarly, the volume integral of interior stress equation (9.3) can 
be expressed as 

f T 

J j B J 1 j k (x.T;I,t)A®J 1 (x.T)dvaT 

o V 

L T 

' 2 I I B7 lj]t [ J1 ,n (l).T;i.t)]M ? ( a )d7 m 4.J 1 ”ar (9.11) 

m=l o V 
m 


in which the time integral is treated analytically as before. 

(B) Spatial Integration 

The nonsingular, spatial integration of volume integrals of equations 
(9.2) and (9.3) are evaluated numerically by applying the Gaussian 
quadrature technique of the transformed integral as 


ill 

j B[x(n),£] Mp(a)dV m = J I j" B[x(a.) ,i]Mp(2.) J(Tj.)dri 1 dii 2 dTi 3 
v m -1 -1 1 


ABC 

= 555 W^Ix,^),!] Hp(r, abc )J(ii abc ) (9.12) 

a=l b=l c=l 


129 



where the Jacobian is defined by 

= J(n 1 ,n 2 .n 3 )dT» 1 dT l2 dn3 
and is given explicitly as 


I 3x. | 

= j 3rjT j = 1,2,3 


For singular volume integrals, the volume cell can be transformed to a 
unit cube and the cube is subdivided into tetrahedra through the field 
point, as shown in figure 9.2. Using a local spherical polar coordinate 
system (r.e.d) with its origin at the field point, the integral of the sub- 
cell can be transformed by the Jacobian as 


2 

dV^ = JdrdbdO = r sina drdddQ 


The integrand involving the kernel is singular of the order l/r 2 and 

therefore the integral is bounded in the transformed domain. However, the 
volume integral is singular of the order l/r 3 and in the 

transformed domain the behavior is approximately of the order l/r . The 
integral, however, is made bounded by excluding a sphere and mapping the 
remainder of the tetrahedra to a unit cube as shown in figure 9.3. The 
integration is computed by applying the Gaussian quadrature to the 
transformed domain. A series of numerical trials with different sizes of 
the spherical exclusion led to the surprising conclusion that it could be 


set to zero for the most accurate three-dimensional analysis. 

The above described volume integration scheme is based on the work of 
Mustoe (1984), Bajernee and Davies (1984), Raveendra (1984) and Banerjee 
and Reveendra (1985). 
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IX. 5 TIME-STEPPING AND ITERATIVE SOLUTION ALGORITHM 
(A) Time-stepping 

In order to obtain the nonlinear transient response at a time T^, the 
time axis is discretized into N equal time intervals, i.e. 

N 

- J nAT (9.13) 

n=l 


where AT is the time step. 

Using equation (9.13), the integral equation (9.2) can be written as 


°ii 4 u i ( i-V - J j [< 3 ii 4t i - F ii 4 u i ldsd ' - J I 


t n-i s 


T V 
i N-l v 


.N-l f N r 

J I tG ij At i ' fij^ildSdr ♦ J j 

T=0 S 


T V 
N-l 


(9.14) 


For the present case, the linear time interpolation scheme described in 
Sec. VIII. 3 .B is used to approximate the time variation of the field 
quantities during a time step because the same scheme can also be used for 
constant time interpolation with averaging. 

Thus, after the usual discretization and integrations (time and 
spatial both), the integral equations (9.14) are transformed into an 
assembled system equation of the form 

o 

[A^HAX 1 *} - [B*](AY N ) - [C*]{Aa N ) 

N 

= ~ \ [[aJ+aJ - 1 ] {AX N_n+1 ) - [bJ+bJ - 1 ] {AY N_n+1 } 

n=2 
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+ [Cj+cJ^HA a N - n+1 l] 


(9.15) 


[A^]{AX N } = (Bj](AY N } + [d]{A<r N } + (R N ) 


(9.16a) 


a, b AX N = Afc bN + Ak° bN 


(9.16b) 


where A and B are the matrices related to the unknown and known 
incremental displacements and tractions; 

C is the matrix related to the initial stresses; 

AX and AY are the vectors of unknown and known incremental 
displacements and tranctions; 

for AX, AY and A<r, superscript denotes time, i.e. AX n = X n - X n 1 ; 
for A, B and C matrices, superscript denotes the time step when 
they are calculated, and the subscript denotes the local 
time node (l or 2); 

R N is the effect of past dynamic history 
A& bN = [Bj] {AY N } + (R N ) ; 

A b = CA^l , and 
Afc° bN - [C*l {A<t N J . 

Similarly, the integral equation for stresses can be written in a 
discretized form as 


(Ao N ) = [A*] (AX N ) + [B*HAY N ) + CC^] CAct N ) + {R aN ) 


or a/ = A c AX N + Ab° N + Ab 0aN 
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where 


# _ 9 


indicates that the matrices are related to the stress equation; 


A ff = [aJ] ; a/ = [B*] {AY N } + {R ffN } ; and 
Afc 0<TN = EC*] {A° N } 


The algorithm described here provides the solution of system equations 

given by equations (9.16) and (9.17). The solution of these system 

equations requires complete knowledge of the initial stress distribution 
°N 

A£ within the yielded region that is induced by the imposition of the 
current increment of boundary loading. This, unfortunately, is not known a 
priori for a particular load increment and therefore an iterative process 
must be employed within each time step. 

This incremental algorithm can be described as follows: 

(i) Obtain the transient elastic solution for an arbitrary increment of 
boundary loading AY N during the time interval T n-1 to T N , as 

- w* 


Ah N - A^Ax' 3 + Ah < ’ s 


where N is the time step number. 

If the material has not yielded yet, accumulate X-vectors, i.e. 
X N = x N_1 + a* n . 


(ii) If the material was yielded before go to step (vi) . 

(iii) Check whether any node has yielded during the current time step. If 
the material has not yielded yet, accumulate stress and strain, and go 
back to step (i). 
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(iv) Calculate the value of ct q , equivalent stress by using ct T = + 

Act W as the stress chnages and compile a list of yielded nodes. For 
elastic nodes accumulate the stress and strain, i.e., ct N = ct T and 

= e. N_1 + lD e ] -1 AcL N . Calculate the correct stress at the elasto- 
plastic nodes by using the elastoplastic stress-strain relations A» e P 
= Q ep Ae. and using the elastic strain increments as a first 
approximation. Modify the stress history for yielded cells ct N = ct^ - 1 
+ Aff ep . Calculate initial stress Act 0 = « T - ct N . 

(v) Assume Ab^ N = 0 and Ak ffN = 0 and using the generated initial 

stress Act 0 calculate a new A& N by using equation (9.16b) and Act N 

by using equation (9.17b). Calculate the equivalent stresses by using 
T N N 

the history ct = ct + Act and compile a list of yielded nodes. For 

elastic nodes, accumulate the stress ct = ct x and strain. For the 

elastoplastic nodes calculate the currect stress Aw ep = Q ep Ae . The 

initial stresses generated are Aa° = Act N - Aw ep . Modify the stress 

history for the yielded nodes ct N = ct N + Aa ep . Accumulate A2£ N and Act N 

o o 

(i.e. A£ N = a^ N + AX^ and Aa^ = Act^ + Act 0 ), so that they can be used 
in the next time step for past convolution. 

(vi) Check if the initial stresses Act 0 are less than the acceptable norm 
and if so go to step (i) and if not go back to step (v). If the 
number of iteration exceeds, say, 50 then it is reasonable to assume 
that collapse has occurred. 

IX. 6. EXAMPLE OF APPLICATION 

In order to demonstrate the accuracy and applicability of the proposed 
nonlinear transient dynamic analysis algorithm, a presentative problem is 
analyzed. English units are used with foot (ft) for length, pound (lbf) 
for force, and seconds (s) for time. 
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(a) Bar subjected to a step end load . 

A bar with circular cross-section is held along its sides by 

lubricated rollers and is fixed at one end. The free end is subjected to a 

suddenly applied and maintained uniform compression t = -333 which 

z 

exceeds the yield stress of the bar (i.e. yield stress of the bar is Y = 
300). in this example, the bar has dimensions and material properties 
identical to that of example VIII.6.b(iii). The discretization of the bar 
is similar to the one shown in figure 8.8 except, in the present example, 
the full cross-section of the bar is modeled instead of one-quarter of it. 
The volume of the bar is discretized by using five 20-noded, volume cells 
of equal dimensions. A bilinear stress-strain relation as shewn in figure 
9.4, is assumed to describe the rod's material property. The time step 
used for this example is AT = 0.004473. In figure 9.4, the elasto-plastic 
response of the bar at time T = 0.8 T g (where T e = c^T/L , i.e. the time 
taken by the compression wave to reach the fixed end of the bar) is plotted 
against the one-dimesnional analytical solution (Ref. Garnet and Armen, 
1975). In this, the normal stress a ^ are normalized by the elastic 
modulus and the distance along the bar is normalized by the length of the 
bar. The numerical results are in reasonable agreement with the analytical 
solution except for the sharp jumps in the stress which are diffused by the 
numerical analysis. The major differences in the results between the two 
solutions can be attributed to the three-dimensional nature of the present 
example. As the bar is on lubricated rollers, in addition to longitudinal 
stress, lateral stresses also exist in the bar. Simple one-dimensional 
theory considers longitudinal stress only and thus, the difference between 
the two solutions. 


135 



IX. 7 COMCLDDIttn REMARKS 


A direct boundary element formulation and its numerical implementation 
for nonlinear transient dynamic analysis of three-dimensional isotropic 
homogeneous or piecewise homogeneous solid has been presented. Due to the 
lack of available solutions for three-dimensional nonlinear transient 
dynamic problems, it was found impossible to compare results for a real 
three-dimensional problem. However, the present algorithm is found to 
produce very accurate results for three-dimensional static nonlinear 
problems by using large time steps, (i.e. when the loading is done slowly). 
Similarly, when a large value of yield stress is selected, the incremental 
nonlinear transient algorithm is found to produce results identical to that 
produced by the linear transient algorithm. This new formulation provides 
a numerical tool for solving three-dimensional transient problems involving 

material nonlinearity which are now impossible to solve by any other 

» 

method. 
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CHAPTER X 





X.l GENERAL CONCLUSIONS 

A complete and general numerical implementation of the direct boundary 
element method applicable to free- vibration, periodic vibration, and linear 
as well as nonlinear transient dynamic problems has been presented. The 
developed methodology is applicable to problems involving two or three- 
dimensional, isotropic, piecewise-homogeneous solids of arbitrary shape. 
Since all of the proposed analyses are based on the boundary element method 
(BEM), they have all the advantages of the BEM over the Finite Element and 
Finite Difference methods such as, discretization of only the boundary of 
the domain of interest rather than the whole domain, ability to solve 
problems with high stress concentrations, accuracy and the ease of solution 
in infinite and semi-infinite mediums. 

The real-variable BEM formulation presented in this dissertation 
provides a numerical tool for free-vibration analysis of solids with 
complex geometries. This method has been compared with MARC-HOST Finite 
element analysis and was found to yield essentially similar results for a 
cantilever beam problem. Thus, the proposed method is a viable alternative 
to algorithms based on Finite element schemes. In addition, it needs only 
the boundary discretization of the problem rather than the whole domain. 

The advanced implementation of the BEM for steady-state dynamic 
analysis of two and three-dimensional, visco-elastic solids, presented in 
chapters IV and VI, are one of the most general numerical implementation 
presently available. By comparing the results obtained by the present 
impl orientation with those by other methods, the accuracy and stability of 
the present method is established. For half-space problems, the proposed 
methodology is a better alternative to the conventional finite element 
method. For half -space problems Finite element presents two restraints: 
(i) the model must be bounded at the bottom by a rigid bedrock, and (ii) 
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the soil away from the vicinity of the foundation is represented by 
parallel layers unbounded on the horizontal direction. These two 
conditions are not always close to reality whereas, in BEM, the fundamental 
solution satisfies the radiation condition at infinity and therefore no 
bounding surfaces are needed and only a small number of elements are 
necessary to model the problem. 

The transformed-domain boundary element formulation presented in 
chapter VII is capable of providing accurate solutions to transient elasto- 
dynamic problems. The accuracy and stability of the present implementation 
are established by comparing the results obtained against the available 
solutions from Finite element. Finite Difference and Time-domain Boundary 
element methods. However, the transformed domain formulation suffers from 
the following defects. 

(i) The transform solution is essentially a superposition of a series 
of steady-state solutions and is therefore applicable only to linear 
elasto-dynamic problems. For nonlinear problems, the solution must be 
obtained in the real time domain. 

(ii) Since the Laplace/Fourier transform casts the entire problem in 
the complex domain, the computer time and storage requirements are 
considerably increased. 

The time-domain boundary element formulation for linear and nonlinear 
transient dynamics presented in chapters VIII and IX eliminate the above 
mentioned problems. The proposed time domain methodology, in conjunction 
with the direct step-by-step integration, provides the transient response 
directly and thus it has been extended for nonlinear problems by using an 
iterative algorithm. Using this method, the transient phenomena during 
early response times, preceding the harmonic steady-state motion, can be 
captured while frequency domain methods are incapable of detecting than at 
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all. In addition, approximations related to the value of Poisson's ratio 
and to the number of modal shapes required in frequency synthesis are 
eliminated. 

The versatility of the proposed time-domain methodology is evident in 
view of the results presented in this dissertation for various three- 
dimensional transient problems. Due to its general character, it can be 
used for solving more sophisticated problems. This algorithm is an 
unconditionally stable implicit time marching scheme and is capable of 
producing accurate results. However, for better accuracy, it is 
recommended that the time step should remain snaller than L/Cj ; where L 
being the smallest distance measured along the surface between two corner 
nodes of an element and c 1 being the propagation velocity of pressure 
wave. 

By taking the material nonlinearity into account, the proposed 
methodology for time-domain nonlinear transient analysis has the potential 
to provide a numerical tool for solving soil-foundation problems in a more 
realistic manner which cannot be accomplished by using the available 
transform domain algorithms. 

X.2 RECOMMENDATIONS 

In order to facilitate future research based on the findings of the 
present work the following are recommended: 

1. The stability of the time-domain transient dynamic algorithm has 
been established for simple problems by analyzing the problem of 
radial expansion of a cavity in an infinite space for different time 
steps and meshes. However, to insure the stability and convergence of 
this algorithm for more sophisticated problems, further investigation 
by using a complex problem is recommended. 
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2. As mentioned earlier, the transformed domain dynamic analysis 
yields erroneous results when the forcing frequency happens to be one 
of the natural frequencies (or fictitious eigenfrequencies in the case 
of exterior problems) of the structure under consideration. To 
eliminate this problem, a computationally feasible modification of the 
transformed domain algorithm is needed. 

3. In the present work, for nonlinear dynamic analysis, the Von 
Mises constitutive relations are used to model the material behavior. 
However, for materials like soils, a more realistic material model 
needs to be included to model the nonlinear material behavior during 
dynamic loadings and unloadings. Moreover, only a simple test problem 
has been solved in the present work. However, for solving realistic 
engineering problems further work is needed. 

4. • The problems of soil-structure interaction during an earthquake 
excitation is of considerable importance to civil engineers. This 
problem can be tackled in a deterministic way by modifying the present 
formulations. For this purpose, extension of the present algorithms 
to solve the general wave scattering problems by including the 
effects of incident waves in the formulation is recommended. 

5. The proposed time-domain transient formulation involving 
convolutions provides accurate results, but it is computationally 
expensive. However, for certain class of problems such as those 
related to structural dynamics, an approximate and computationally 
inexpensive boundary element formulation can be developed. This can 
be achieved by extending the method proposed for free-vibration 

s 

analysis to linear and nonlinear transient dynamic analysis of solids. 
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6 . In practice, inhomogeneity and anisotropy are present in most 
engineering problems. Whilst the inhomogneity can be handled by 
substructuring, it is of extreme importance to develop appropriate 
fundamental solutions for dynamic analysis of problems involving 
anisotropy . 

7. Some of the dynamic problems such as non-destructive testing of 
materials involve material nonlinearity as well as geometric 
nonlinearity. Therefore, extension of the present nonlinear transient 
dynamic formulation to include geometric nonlinearity is desirable. 
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Two-dimensional boundary elements 
Figure 4.1 


156 


157 



Boundary Element discretization of a half-space Problem 


Figure 4.2 
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Figure 4.4 
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Figure 4.5 
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REAL PART OF CONTACT STRESS FOR VERTICAL 
VIBRATION OFA RIGID STRIP FOOTING 

Figure 4.6 
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IMAGINARY" PART OF CONTACT STRESS FOR 
VERTICAL VIBRATION OF A RIGID STRIP FOOTING 
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IMAGINARY PART OF CONTACT STRESS FOR 
HORIZONTAL VIBRATION OF A RIGID STRIP FOOTING 


Figure 4.9 





REAL PART OF CONTACT STRESS FOR 
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Figure 4.10 
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IMAGINARY PART OF CONTACT STRESS FOR 
ROCKING OF A RIGID STRIP FOOTING 


Figure 4.11 
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DISCRETIZATION OF A MACHINE FOUNDATION 
ON AN ELASTIC HALF SPACE 


Figure 4.12 
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IMAGINARY PART OF STIFFNESS CO-EFFICIENTS 
FOR A MACHINE FOUNDATION 

Figure 4.14 
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REAL PART OF STRESSES FOR VERTICAL VIBRATION 
OF A MACHINE FOUNDATION 
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VIBRATION OF A MACHINE FOUNDATION 
Figure 4.16 
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REAL PART OF STRESSES FOR ROCKING OF 
A MACHINE FOUNDATION 


Figure 4.17 
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IMAGINARY PART OF STRESSES FOR ROCKING OF 
A MACHINE FOUNDATION 
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Boundary element discretization of a cantilever beam 

Figure 5.3 
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Discretizations of a shear wall 
Figure 5.4 
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Figure 5.5 




Three-dimensional nonplanar surface patch 
Figure 6.1 
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(b) Eight-noded quadtrilateral 


Three-dimensional surface elements 
Figure 6.2 
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Infinite element 
Figure 6.3 
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Singular integration process for a quadrilateral element 
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Cantilever subjected to harmonic end shear 
Figure 6.6 
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Figure 6.7 
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Time history of displacement u 2 at the internal point E( 0,-60') 
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AXIAL DISPLACEMENT AT THE FREE ENO 

Transient Analysis of a Cantilever Subjected to a Harmonic 
Axial Loading 
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Time-marching process 
Figure 8.1 
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Normalized radial displacements of the cavity surface 
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Deviator ic stresses at the cavity surface for suddenly 
applied and maintained pressure 
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Hoop stresses at the cavity surface for a cavity 
engulfed by a pressure wave 
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Radial scattered displacements for a cavity engulfed by a pressure wave 
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Boundary element discretized for a point load 
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Three-dimensional volume cell 
Figure 9.1 
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Geometrical mapping of subcell (excluding spherical segment) 
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Figure 9.3 
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APPENDIX A1 


BOUNDARY KERNELS FOR TWO-DIMENSIONAL STEADY-STATE DYNAMICS 


The tensors and Fy in the transformed domain are of the form: 


Gy <!.!..> - JJJ tASy - Bt^ty] 


(Al.l) 


= 2n [p(5 ij an + r ,j n i ) 


+ Q(r -n- - 2r 
• j 


3r 3r , 

i r -i + ^ i r i + Sr i r 

l ,3 *3 *3 


(A1.2) 


where 

A = K q (z 2 ) + (l/z 2 ) [K 2 (z 2 ) - (c 2 /c 1 )K 1 (z 1 )1 
B » K 2 (z 2 ) - (c 2 /c 1 ) 2 K 2 (z 1 ) 

P = 3A/ar - B/r 
Q = - 2B/r 
R = - 23B/3r 

S = {(c 1 /c 2 ) 2 - 2} (3A/3r - 3B/3r - B/r) 

n = normal vector (Al.3) 


where z 1 = imt/c 1 and z 2 = iwr/c 2 ; and 

K q , Kj and K 2 are the Modified Bessel functions of second kind, 
having the following recursive properties: 


A-l 


K 2 (z) = K 0 (z) + (2/z)K 1 (z> 


K'(z) = - K. (z) 

0 1 

K:(z) = - K. (z) - K- (z) /z (A1 .4) 

1 o 1 

where the bar denotes the differentiation w.r.t. z . 

Using the recursive formulas (A1.4) the tensors and can be 

expressed in terms of Modified Bessel functions of the second kind of 
orders zero and one. These functions are given below along with their 
expansions for small and large arguments. 

Modified Bessel Function of Second Kind 
Zero order : 


K 0 (z) = tln(z/2) + rll 0 (z) + 


z 2 /4 

( 1 !) 


2 + (1 + 1 / 2 ) 


( z 2 / 4 ) 

( 2 !) 


2 

2 


+ (1 + 1/2 + 1 / 3 ) 


(z 2 /4) 3 
( 3 ! ) 2 


(A1.5) 


I o ( z) 


, * z 2 / 4 . ( z 2 / 4) 2 ( z 2 / 4) 3 (z 2 / 4) 4 

(l!) 2 (2! ) Z (3!) Z (4! ) Z 


(A1.6) 


where r - 0.5772156649 


First order : 

v (z 2 /4) m 

K x (z) = (l/z) + ln(z/2)I 1 (z) - z/4 2 tMm+l) + Mm+2)] ^^ 77 

m=o 

(A1.7) 


A-2 



_ , 2,vm 

x i< z > " <*/« 2 iTF(iT2T 

m=o 

where <l) *» — r » and 

n-l 

(n) * - r + ^ m_1 for n 2 2 
m=l 


Snail argument ■ expansion : 

If z -> 0 (i.e. abs(z) < 10 -8 ) 

K q (z) = - ln(z) and K^z) = 1/z 


Large argument expansion: 

If z is large (i.e. abs(z) > 3.5) 


K c (z) 



(l) 2 , (3) 2 _ (3) 2 (5) 2 

8z 2! ( 8z) 2 3!(8z) 3 


+ (3) 2 (5) 2 (7) 2 
4! (8z) 4 


] 


k i (2 > - Jk e ' Z t 1 + fe - 


3x5 + 3x5x21 

2!(8z) 2 3 ! ( 8z ) 3 


However, for abs(z) > loo , K Q (z) = Kj^z) o.o 


(A1.8) 


(A1.9) 


(A1.10) 


(A. 11) 


(A1.12) 


(A1.13) 
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APPENDIX A2 


The tensors G^j and in the transformed domain are of the form: 


(x.i.s) = 4ff|I [A5 i; . - Br ^r 


(A2.1) 


(x.i.s) 


- h [ p(5 ij i * c .j n i> 


+ Q(n-r . - 2r -r • |§) + Rr -r _= + Sr -r .-l 

j ,1 .1 o a n >1 *] ®n o 


'(A2.2) 


where s is the Laplace transform parameter. In addition. 


3c, 3c, 

( + — + l ) 

s 2 r 2 sr 


" Sr/C 2 2 2 

> c, 3cf 3c, -sr/c. 

i < “A + — + 1 ) e 

r c 2 s 2 r 2 sr 


(A2.3) 


C 2 C 

< 4 V“ + 1 > 

s 2 r 2 sr 


" SC/C 2 2 2 „ - Sr/C l 

: 

r c\ sV sr 


(A2.4) 


where e is the exponential function. Furthermore, 


p_M B q__ 2B n _ _ 2 — 

p 3r r ' Q r ’ R 2 3r ' 


_2 9r 3r r 

c 2 


(A2.5) 
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APPENDIX A3 


TMTERTOR STRESS KERNELS FOR STEADY- STATE DYNAMICS 
The interior tensors and F^j k for two-dimensional steady 

state dynamic analysis are of the form: 


3G . 3Gm 3Gi. • 

- « jk * Jj-) 


(A3 .1) 


_ 3F . 3F . . 3F t . 

F ijk < *'*'* ) = l5 jk + * Sif’ 


(A3. 2) 


where 


3G ij 

a ?k 


2jtji [ 5 ij 3r r .k “ 3r r .k r .i r ,j B(r .i r ,;jk + r ,j r ,ik ) ] 


(A3. 3) 


3F. . 


ii l rap ar 

it 1 - - u L« r ik (s tj «5 * r .j n i ! * !><6 y r ,mk n in + r , jk n i’ 


L 3Q * 3r. . 3r. 

+ 37 *■ > k r » i°j - 2r .i r ,j 3i ) + ^.ik^ - 2r ,ik r ,j in 


2r .i r ,jk 3 q " 2r ,i r ,j r ,mk n m + 3r r ,k r ,i r ,j 3n 
+ R * r .ik r .j 3n + r .i r .jk 3n + r . i r » j^mk 1 ^ 


+ 


as 

3r 



+ S r 



(A3. 4) 


where the functions A,B,P,Q,R and S are listed for two and three 
dimensional problems in Appendix A1 and A2 , respectively. 
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The tensors and F^j are of the form: 


Gyd.TSl.O - ^ [(SBy-by) J «<V«>d>. 

i/c. 


+ a^td/c JSfv-r/Cj) 


- <l/c 2 )6(v-r/c 2 )} + (b^/c^^tv-r/cj)] 


3 

a.. = y.yj/r and b i;j = 5^/r 
" x i 


P ii (X.T;i.r) = ^ [-ecjcsay-by) J XS(v-Xr)dX + (da^-Ib..) 


{5(v-r/c 2 ) - (c 2 /c 1 ) 2 8(v-r/c 1 )} + Zra^/CjCS ' (v-r/c 2 ) 


- ( c 2 /c 1 ) 3 5'(v-r/c 1 )}- c^. (l-2c 2 /c 2 ) (SCv-r/Cj) + (r/CjJS'tv-r/Cj)} 


- d^. (5(v-r/c 2 ) + (r/c 2 )S'(v-r/c 2 )] 


where v = T - r 


a ij - yjyjW' 

c ij - yj n i /r3 

d ij = (y.n. + 8 ij y m n m )/r 2 


b ij " c ij + d ij 
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INTERIOR STRESS KERNELS FOR TRANSIENT DYNAMICS 


The tensors G?j k and F®j k are of the form: 

l/c 2 

Gf jk (X.T :i .«) - ij [&£<» ljk -b iJk > J I5(v-U)d>. 

l/c. 


- (I2a ijk -2b i;jk ){8(v-r/c 2 ) - (c^Cj^v-r/c^ } 

- 2ra ijk /c 2 {5'(v-r/c 2 ) - (c 2 /c 1 ) 3 8'(v-r/c 1 )} 

+ c ijk (1_2c 2 /c i ) f5(v-r/c i ) + ( r /c i )8 ' (v-r/c i )} 

+ d i j k (5(v-r/c 2 ) + (r/c 2 )5'(v-r/c 2 )}] 


where v » T - t 

a ijk • Wk ,tS 

c ijk “ 6 jk^ i/ r 

d ijk ■ * 6 ijyk )/c3 

^ijk = c ijk + d ijk 

p 1 ^ C 2 

Fj jk <!£.T i l.T> - - «; [l2c5<35a ljk -5b ijk *c ljk ) j X5(v-kr)dk 


1/c, 


4C 2 (45a ijk- 5b ijk +C ijk ){6(v - r/C 2 ) - 


(c 2 / c i) 


'Sfv-r/Cj) ] 
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(A5.1) 


- 4c 2r (ioa i j k -b i j k H5'(v-r/c 2 ) - (c 2 /c 1 ) 3 6'(v-r/c 1 )} 

+ 2c 2 (i-2c 2 /c 2 )(3d ijk -2e ijk H5(v-r/c 1 ) + (r/c 1 )6 ' (v-r/c 1 ) } 


+ c 2 ( ^ f ij k " 2g i: j k ){S(v-r/c 2 ) + (r/c 2 )5 ' (v-r/c 2 ) } 

- 4a i j k r 2 {6"(v-r/c 2 ) - (c^c^V ' (v-r/c^)} 

+ c 2 (l-2c 2 /c 2 ){2c 2 d ijk /c 2 + (l-2c 2 /c 2 )e ijk )8"(v-r/c 1 ) 

* t2f i:k 6 ' ,<v " c/c 2 > ] 

where v = T - r 

a ijk • ¥jVA /[7 
d ijk - ( yj5'k n i + 
e ijk = 8 jk n i/ c 

£ ijk - ‘y^k + y^j + ymV 5 ikyj* 8 ijyk ))/rS 

9ijk ■ (6 ijV S ik n j l/c3 

b ljk “ d ijk + £ ijk 
c ijk = e ijk + g ijk 
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VOLUME KERNELS FOR TRANSIENT DYNAMICS 
The tensors and jJj^ are as follows: 

l/c 2 

B ilj <x,T 1 i.t) = ^ [-(lSa^-Sb^j) j W(v-lr)dl * l/c|(«a llj -b llj ) 

1/Cj 

{8(v-r/c 2 ) - (c 2 /c 1 ) 2 8(v-r/c 1 )} + ra^j/CjfS' (v-r/c 2 ) - (c 2 /c 1 ) 3 5'(v“r/c 1 )} 


- c il j/2c 2 {8(v-r/c 2 ) + (r/d 2 )8’(v-r/c 2 )}] * 


(A6.1) 


where v = T - t 
a uj ■ 

°uj ■ ‘Vi * 5 ijn )/t3 
b uj - c uj + s uyj /r3 


B? 


iljk 




-U* 


2 (35a iljk" 5b iljk +c iljk ) 


l/c 2 

I 

1/c, 


X8(v-Xr)dX 


- 2 ( 45a n jk _6fcl il jk^Cil jk> t 5( v-r/c 2 ) - (c 2 /c 1 ) 2 8(v-r/c 1 ) 

- 2 r/ c 2 < 10a n ^ 5 ' (v— r/ c 2 > - (c 2 /c 1 ) 3 5'(v-r/c 1 )} 

+ (3f iljk -g iljk ) (1-2c 2 /c l ) f6(v “ r/c l ) + (r/c i )6 ' (v-r/c i )} 
+ (l.5d^jk ~ eiij k ) {8(v-r/c 2 ) + (r/c 2 )8'(v-r/c 2 ) 
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- (2r2 a il j k /c 2 ) (5"(v-r/c 2 ) - (c 2 /c 1 ) 4 5"(v-r/c 1 )} 

+ (r 2 f il j k /c 2 ) (1-2 c 2 /c^) 5" (v-r/Cj) + (r 2 d iljk /2c 2 )6" (v-r/c 2 ) ] 

where v ■ T - x 

a iljk * Wk' 1 ' 

a ujk ■ 'Vi^k * 'l/ft * 6 ik y i y j * 6 ikyiyj ,/r5 

3 

e il;jk = (s ij 5 lk + 8 ik 5 lj )/r 
£ ujk • 

9 iljk * 8 jk 8 U /c3 

b iljk = d iljk + f iljk + 5 il y j y k /r5 
c iljk = e iljk + g iljk 

F qx. J - D : 


J iljk 


15(1- v) 


[(7-5v)8.jJ 


kl 


+ <l-Sv) s i iS j k] 


(A6.3) 
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APPENDIX B 


PROPAGATION OF WAVEFRONTS AS SURFACES OF DISCONTINUITY 

When a body is disturbed from a quiescent state by excitation at a 
portion of the boundary or within a restricted domain inside the body, 
neighboring domains are soon set in motion and put into states of 
deformation. The moving surface which separates the disturbed from the 
undisturbed part of the body is called the wavefront. At wavefronts, the 
field quantities and/or their derivatives may be discontinuous. However, 
if the material remains coherent and does not fracture, the displacements 
will certainly be continuous in both space and time. In many situations, 
involving very sudden loadings, the particle velocities and/or stresses 
will have sudden variations (discontinuities) at the wavefront over a very 
small interval of space and time. These variations at the wave fronts can 
be quite closely approximated by finite jumps based on the basic techniques 
developed towards the end of last century for the study of propagating 
surfaces of discontinuity in continuum mechanics. 

Love (1904) sets down the following basic kinematical and dynamical 
conditions that must hold at a propagating surface of discontinuity in an 
elastic solid. 

Kinematic conditions : 

Consider a surface of discontinuity S , propagating in an unbounded 
medium. The situation is shown in figure B.l, for a fixed instant of time. 
It is assumed that S propagates into region (2), leaving a region (l) 
behind it, and moves normal to itself with velocity c , i.e., each point 
P(x) of S propagates with velocity c , along the outward unit normal 
vector q to S at that point. If one supposes the components of <*„ 
are discontinuous across S , the jumps will be denoted by the standard 
bracket notation. 
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s 



where the subscript 1 denotes the value of the field variable on S when 
S is approached through region (l), and the subscript 2 is employed to 
denote the value when S is approached through region (2). 

As mentioned earlier, since the material should maintain its integrity 
at the wavefront, the jump in the displacement components at S is zero, 
i.e. 

lu-^3 - 0 (B.2) 

Moreover, if the strains and velocities at a wavefront are discontinuous 
(i.e. shock waves), the finite jumps in them must satisfy the following 
kinematic relations. 

• 

Eu. ] + cEu. .]n. = 0 
i J j 

Eu^] + c[3u i /3n] = 0 
• 

Eu-ln- + cEu- = 0 (B.3) 

i j it j 
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However, if the first derivatives of the displacements across S are 
continuous but the second derivatives are discontinuous (acceleration 
waves), then the following kinematic relation has to be satisfied at the 
wavefront . 


'“i.jk 1 - 'v’j n k - 0 

t(u i ) j] + cq^Rj = 0 

Cu i 3 - c 2 q i = 0 (B.4) 

where q^ is an unknown function. 

Dynamical, conditions : 

The dynamical conditions, which has to be satisfied at the moving 
surface of discontinuity S , are determined by considering the momentum 
changes of a thin slice of the medium adjacent to S and the corresponding 
impulse-momentum equation. It has the form 

t^ijlnj + pc[ Ui ] = o (B.5) 

For acceleration waves, the jumps in the second derivatives of u should 
satisfy the linear momentum equation, i.e. 


X[u . ] + |i [u. .. + u- ••] = p[u.] 
m,im J ^ i,]] j.rj i 


(B.6) 


The fundamental singular solution of transient elastodynamics, for the 
displacements generated by a suddenly applied concentrated load at a point 
of the unbounded elastic medium was first developed by Stokes (1899). 
Love (1903) performed an extensive study of Stokes' solution for initial 
value problem with arbitrary initial values, and related wavefront 
discontinuities. He pointed out that Stokes' formula yields correct 
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results only when the input field quantities are continuous at the 
wavefront. He also found that the Stokes' formula satisfies the necessary 
continuity conditions on the displacements (eq. B.2), the kinematical 
conditions on the velocities and strains (eq. B.3) and the dynamical 
conditions on the stresses and velocities (eq. B.5), provided the input 
function is continuous. Thus, if the input excitation is a step loading, 
it has to be modeled as a ramp loading in the first time step. Also, a 
very small time step cannot be used for this purpose, because it will 
result in non-vanishing dilation and rotations at the wavefronts. 
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APPENDIX Cl 


TSOPARAMETRTC BOUNDARY ET.EMEOTS FOR 2-D PROBLEMS 
Both the three noded quadratic and two noded linear elements were 
depicted in figure 4.1. The shape functions for the three noded quadratic 
elements are: 

N^i,) = 2(1, - 1/2) (i, - l) 

Nj (n ) = -4 tj (i, - 1) 

N 3 (t,) = 2i,(i, - 1/2) (Cl.l) 

where i, is the intrinsic coordinate (0 <. r, <. l). 

The shape functions for two noded linear element are: 

N.^1,) =>1-t, 

N 2 (i,) = i, (Cl. 2) 

The normal unit vectors along the positive x and y axes are defined as: 
n 3 = Oy/3i,)/|j(i,) I 

n 2 = (-3x/3i,) |j(n) I (C.1.3) 

where |j(i,)l is the magnitude of the determinant of the Jacobian matrix 
(Ref. Sec. IV. 4) . 
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APPENDIX C2 


Both the six node triangular and the eight node rectangular surface 
elements were depicted in figure 6.2. It is worth noticing that these 
intrinsically planar elements becomes curved in three-dimensional space. 
Hie shape functions for the six node triangle are: 


n (2 ti - 1) if a = 1.2.3 
a a 


4 Vr 


if a = 4.5.6 with p = a - 3 
and y = a - 2 


(C2.1) 


where and n 2 are two linearly independent coordinates and ti 3 = 
The shape functions for the eight node rectangle are: 


0.25(l+? o )(l+n o )(5 o+T , o -l) if a - 1,3. 5.7 


N a = 0.50(1-5*) (l-ii 0 ) 

0.50(i+5 o )(i-n 2 ) 


if a = 2,6 
if a = 4.8 


(C2.2) 


where 5 Q = 5 5 o and n ii o , with 5 and n being the two linearly 

independent coordinates and (5 ) the coordinates of node a . 

cl a 

Two base vectors along the intrinsic coordinates 5,n (or can 

be defined as 


9x . 3v t 9z , 

£ 1 “ 35 ^ ^ + 35 i + 95 ^ ^ 


9x , . dy , . 9z , , 
£ 2 = dr, i + dn a + ^ d n k 


(C2.3) 


where i, j, and is are unit vectors along the x, y, and z coordinates. 
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respectively. Their cross product 


% = x £ 2 (C2.4) 

is a vector normal to the surface of the element and its magnitude is equal 
to the value of the determinant of the Jacobian matrix. 
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APPENDIX D1 


ANAT.YTTCAL TEMPORAL INTEGRATION OF THE TRANSIENT 
DYNAMIC KERNEL FOR CONSTANT TIME INTERPOLATION 


For constant time interpolation the field variables are expressed as 


N 


f(x.t) = ^ f n <x)<a n <T) 

n=i 


(Dl.l) 


where 

d n (T) =* [H(r - (n-l)AT) - H(t - nAT)l ; 

f n (x) represents the spatial variation of the field variable 
f(x,r) at time T n ( = nAT) ; 

N is total number of time steps; and 
H is heaviside function. 


Each of the transient dynamic kernels listed in Appendices A4, A5 and A 6 
has one or more of the following time functions: 


(1) 

8(T - x - r/c) 

(D1.2) 


1/c 2 


(2) 

J X6(T - -c - r/c)dx 

(D1.3) 


Hc 1 


(3) 

• 

8(T - x - r/c) 

(D1.4) 

(4) 

5(T - r - r/c) 

(D1.5) 

where c 

is either pressure wave velocity Cj 

or shear wave velocity c. 


and 8 is the delta function. 

Using equation (Dl.l), the time integrals related to the above time 
functions can be expressed as 


r T 

J g(T - x - r/c)f (x.t)dT 
o 
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(D1.6) 


N nAT 

= ^ [f n (x) J g(T - t - r/c)6 n (-c)dt] 
n=l (n-l)AT 

The time integrals on the right hand side of equation (Dl.6) are 
evaluated analytically as follows. 


Time f u nction i ; 

nAT 

| 5(T - r - r/c)d n (r)dT = <6 n (T - r/c) 

(n-l)AT 


Time function 2 : 


nAT l/c„ 


j J X5(T - r - Xr)dXd n (r)dt 
(n-l)AT l/c. 


l/c 2 nAT 

J J X8(T - x - Xr)d n (r)dtdX 

1/ (n-DAT 


r /<=2 r 2 i 1/c 2 

J Xi„(T - Xr)dX - [a 2 /2)d„<T - It) J x 2 

1/C, 


(D1.7) 


(D1.8) 


An important characteristic of the transient dynamic kernels is the 
time translation property (Ref. Chapter VIII and Appendix A4). Because of 
this characteristic, at each time step only the effects of the dynamic 
history of the first time interval on the current time node needs to be 
evaluated; i.e. at each time step the analytical time integrations has to 
be done only for n = l . Thus, equations (Di.7) and (Di.8) reduce to 

at 

J 8(T - t - r/Od^Odr = ^(T - r/c) = H(T - r/c) - H(T - AT - r/c) 

0 (D1.9) 
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r 4 T r /C2 r 2 -t 1/c 2 

J J XS(T - t - r/cJdXd^TJdr = [(X 2 /2)«5 1 (T - Xr)] 1/c 2 

. « 1 


o 1/c 


.l/c 


= [(X 2 /2)H(T - Xr) ] 1/c 2 - [<X 2 /2)H(T - AT - xr)] 


1/C 2 

l/c. 


(D1.10) 


where 


[(X 2 /2)H(T - Xr)] 


1/c 2 

l/c. 


= 0 


if T < r/Cj . 


2r 2 2c 2 


if T > r/Cj and T < r/c 2 


_1 1 _ 

2c 2 


if T > r/c„ 


(Dl.ll) 


The second term on the right hand side of equation (Dl.9) can be 
obtained in a similar manner by replacing 'T' by 'T-AT' in equation 
(Dl.ll). 

Time function 3: 

The time integrals involving time function (3) are approximated by 
using a backward finite difference scheme, i.e. 

nAT 

J 8(T - x - r/c)f < 2 £,r)dt 
(n-l)AT 


nAT 


J 8 (T - x - r/c)f (x,T)dT 
(n-l)AT 


£ n ( & } ~ f n-l ( * } 


AT 


d R (T - r/c) 


(D1.12) 
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APPENDIX D2 


ANALYTICAL TEMPORAL INTEGRATION OP THE TRANSIENT DYNAMIC KERNELS 
FOR LINEAR TIME INTERPOLATION 


Assuming the field variables to vary linearly during a time step, a 
field variable f(x,r) can be expressed as: 


N 


f(X.T) = } CM 1 f n _ 1 ( X ) +^f n ( x> 1 
n=»l 


(D2.1) 


where and Mj are the temporal shap functions, and are of the form: 




m 2 = 


x - T , 

-TT- V*> 


(D2.2) 


As mentioned in Appendix Dl and Chapter VIII (Section 3), at each time 
step only the effects of the dynamic history of the first time interval on 
the current time node needs to be calculated. Therefore, for n = 1 (i.e. 

T n = AT ’ Vl = 0) 


f(X.-c) = (1 - ^)d 1 (r)f 0 ( X ) + 


(D2.3) 


The analytical time integrations related to ^(t) are same as those 
described in Appendix Dl and the time integrations related to v d 1 (-c) are 
presented as follows. 
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Time-function 1 (eg. Dl.2) : 


-AT 

J 5(T - r - r/cJ-na^TMt = (T - r/Od^T - r/c) (D2.4) 

o 


Time-function 2 (eg. Di.3), : 


l/c„ 


n 

0 l/c. 


X8(T - t - XrJdXrdj^rJdr 


1/c- 


=| X(T - XrJd-^T - Xr)dX 
l/c. 


l/c 2 l/c 2 

= J TXd x (T - Xr)dX - J rX 2 d 2 (T - Xr)dX 
l/c 1 , l/Cj^ 

- T[(x 2 /2)d 1 (T - X/r) ] 1/c 2 “ r[(X 3 /3)d 1 (T - Xr)] 1/( , 2 (D2.5) 


where the term 


[(X 2 /2)d 1 (T - X/r)] 



is evaluated in the same way as 


described in Appendix Dl, and the second term is evaluated as follows. 


[(x 3 /3)d 1 (T - x/r)] 


l/c 

l/c 


2 

1 


l/c 


« [ (X 3 /3)H(T - 5c/r)] 1/c 2 - [(X 3 /3)H(T - AT - xr)] 


l/c, 

4 

l/c. 


(D2.6) 


where 
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[(X 3 /3)H(T - Xr)]j 


if T < r/c. 


T 1 


= — - r if T > r/c.. and T < r/c- 

. 3 - 3 1 * 


3r 3 3c 3 


_l i_ 

3c 3 3c 3 


if T > r/c. 


(D2.7) 


The second term on the right hand side of equation (D2.6) can be 
obtained in a similar manner. 


J 8(T - x - r/c)f (x,f)dT = 


f n (x) " f n-l (x) 


d n (T - r/c) (D2.8) 


(n-l)AT 


The temporal integrations involving the time function 4 (i.e. 5(T - x 
- r/c)) are approximated by using a backward finite different scheme as 
follows: 

nAT 

J 6(T - x - r/c)£(&,x)dx 
(n-l)AT 

nAT 

= j* 6(T - x - r/c)f (x.T)d-c 
(n-l)AT 


f n <x> - 2f n -i ( X> + f n-2 ( * } 


d n (T - r/c) 


(D2.9) 
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